— Z 111  UJ  0 c 

SfZifSOVffl 


E 


0SCCPI  REPQ£TA20 


UNIVERSITY  OF  SOUTHERN  CALIFORNIA 


SEMIANNUAL  TECHNICAL  REPORT 

William  K.  Pratt 
Project  Director 


Covering  Research  Activity  During  the  Period 
1 March  1975  through  31  August  1975 


30  September  1975 


Vy)  Image  Processing  Institute 
University  of  Southern  California 
^ University  Park 

Los  Angeles,  California  90007 

Sponsored  by 

Advanced  Research  Projects  Agency 
Contract  No,  F0606-72-C-0008 
ARPA  Order  No.  1706 


IPI 


The  views  and  conclusions  in  this  doounsnt  are  those  of  the  astbors 
and  should  not  be  interpreted  as  necessarily  representing  the  official 
policies,  either  expressed  or  iaplied,  of  the  Advanced  Bess arch 
Projects  Agency  or  the  O.S.  Governsent. 


I 

If 

T 


OSCIPI  RESORT  620 


SEMIANNUAL  TECHNICAL  REPORT 
Covering  Research  Activity  During  the  Period 
1 March  1975  through  31  Augu&t  1975 


Hiliiaa  K.  Pratt 
Project  Director 
(213)  746-2694 


Image  Processing  Institute 
University  of  Southern  California 
Vniversity  Park 

Los  Angeles,  California  90007 


30  Septeaber  1975 


\ 


^ \ 


This  research  was  supported  by  the  Advanced  Research 
Projects  Agency  of  the  Departaent  of  Defense  and  was 
aonitoced  by  the  Air  Force  Avionics  Laboratory, 
Wright-Pat terson  Air  Force  Base  under  Contract  No* 

P08606-7  2-C-0CC8,  ARPA  Order  No.  1706 


UNCLASSIFIED 

Security  Classification  _____  

DOCUMENT  CONTROL  DATA  - R & D 

(Security  cf»»  Atication  of  tftlt.  body  of  Hbxtrnct  iinti  in  Je  ting  .mnotHtinn  mu  at  h»  ontvrnd  when  thr  ov-trail  report  is  r lusaitiaifj 

I.  ORIGINATING  ACTIVITY  (Corporal*  author)  ^r-  *«.  RF.PORT  SECURITY  CLASSIFICATION 

Image  Processing  Institute  UNCLASSIFIED 

University  of  Southern  California,  University  Park  croup  ' “ ” 

Los  Angeles,  California  90007 


P O R T TITLE 

v^mAGE  PROCE 

f i ^ 


PROCESSING  RESEARCH „ 

* P 


4.  DESCRIPTIVE  NOTES  (Typ»  of  report  And  Inctunv  dmfa) 

Technical  Semiannual,  1 March  1975  through  31  Augug; 


xl 


William  K. /Pratt 


^)iy86O6-72-C-Q(098 


JNARPA  Order  17^6 


Ccknii?<*A.  {tenYiair^  u<M  re^t,  J H\oa.»- 

~ i 31  f7  ^ 


P^RT  

, A 

TOTAL  MO-  STPAjEJ 

seirrai»^75 

_ 4 

Sa 

191 

<GINATOr»  REPORT  NUMBgR(S| 

^USCIPI  6Zff 


9b.  OTHER  REPORT  NO(9>  (Any  other  numbara  that  mmy  bm  maalgnad 
this  report) 


»0.  DISTRIBUTION  STATEMENT 

Approved  for  release;  distribution  unlimited 


II.  SUPPLEMENTARY  NOTES 


I I J.  ABSTRACT 


12.  SPONSORING  Ml  LI  TARY  ACTIVITY 

Advanced  Research  Projects  Agency 

1400  Wiison  Boulevard 

Arlington,  Virginia  22209  


This  technical  report  summarizes  the  image  processing  research  activities  per- 
formed by  the  University  of  Southern  California  during  the  period  of  1 March  1975 
to  31  August  197  5 under  Contract  No.  F08606-72-C-0008  with  the  Advanced  Research 
■N^ojects  Agency,  Information  Processing  Techniques  Office. 

The  research  program,  .entitied^  ^mage  Processing  Research.^has  as  its  pri- 
mary  purpose  the  analysis  and  development  of  techniques  and  systems  for  efficiently 
generating,  processing,  transmitting,  and  displaying  visual  images  and  two  dimen- 
sional data  arrays.  Research  is  oriented  toward  digital  processing  and  transmission 
systems.  Five  task  areas  are  reported  on:  (1)  Imarge  Coding  Projects:  the  investiga- 
tion of  digital  bandwidth  reduction  coding  methods;  (2)  Image  Restoration  and  Enhance- 
ment Projects:  the  improvement  of  image  fidelity  and  presentation  format;  (3)  Image 
Data  Extraction  Projects:  the  recognition  of  objects  within  pictures  and  quantitative 
measurement  of  image  features;  (4)  Image  Analysis  Projects:  the  development  of 
quantitative  measures  of  image  quality  and  analytic  representation;  (5)  Image  Proc- 
essing Systems  Projects:  the  development  of  image  processing  hardware  and  software 
support  systems.  
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ABSTRACT 


This  technical  report  summarizes  the  image  processing  research 
activities  perforaed  by  the  University  of  Southern  California  during 
the  period  of  1 March  1975  to  31  August  1975  under  Contract  No. 
FC86C6-72-C-0CC6  with  the  Advanced  Research  Projects  Agency, 
Information  Processing  Techniques  Office. 

The  research  program,  entitled,  "Image  Processing  Research, » has 
as  its  primary  purpose  the  analysis  and  development  of  techniques  and 
systems  for  efficiently  generating,  processing,  transmitting,  and 
displaying  visual  images  and  two  dimensional  data  arrays.  Research  is 
oriented  toward  digital  processing  and  transmission  systems.  Five 
task  areas  are  reported  on:  (1)  Image  Coding  Projects:  the 
investigation  of  digital  bandwidth  reduction  coding  methods;  (2)  Image 
Restoration  and  Enhancement  Projects:  the  improvement  of  image 
fidelity  and  presentation  format;  (3)  Image  Data  Extraction  Projects: 
the  recognition  of  objects  within  pictures  and  quantitative 
measurement  of  image  features;  (4)  Image  Analysis  Projects:  the 
development  of  quantitative  measures  of  image  quality  and  analytic 
representation;  (5)  Image  Processing  Systems  Projects:  the  development 
of  image  processing  hardware  and  software  support  systems. 
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1.  Research  Project  Overview 
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This  report  describes  the  progress  and  results  of  the  University  j 

1 

of  Southern  California  image  processing  research  study  for  the  period 

j 

of  1 March  1975  to  31  August  1975.  The  image  processing  research 

study  has  been  subdivided  into  five  projects: 

Image  Coding  Projects  . 

Tsage  Restoration  and  Enhancement  Projects 

Image  Data  Extraction  Projects 

Image  Analysis  Projects 

Image  Processing  Systems  Projects 

In  image  coding  the  orientation  of  the  research  is  toward  the 
development  of  digital  image  coding  systems  that  represent  monochrome 
and  color  images  with  a minimal  number  of  code  bits.  Image 

j 

restoration  is  the  task  of  improving  the  fidelity  of  an  image  in  the 
sense  of  compensating  for  image  degradation.  In  image  enhancement, 
picture  manipulation  processes  are  performed  to  provide  a more 
subjectively  pleasing  image,  or  to  convert  the  image  to  a form  more 
amenable  to  human  or  machine  analysis.  The  objectives  of  the  image  ] 

data  extraction  projects  are  the  registration  of  images,  detection  of  j 

objects  within  pictures  and  measurements  of  image  features.  The  image 

l 

analysis  projects  comprise  the  background  research  effort  into  the  3 

i 

basic  structure  of  images  in  order  to  develop  meaningful  quantitative 
characterizations  of  an  image.  Finally,  the  image  processing  systems 
projects  include  research  on  image  processing  computer  languages  and 
the  development  of  experimental  equipment  for  the  sensing,  processing, 
and  display  of  images. 


The  next  section  of  this  report  summarizes  some  of  the  research 
project  activities  during  the  past  six  months,  section  2 is  a list  of 
publications  by  project  members.  Sections  3 to  6 describe  the 
research  effort  on  the  projects  listed  above  during  the  reporting 
period. 


2.  Publications 


The  following  is  a list  of  papers,  articles,  and  reports  of 
research  staff  eeabers  of  the  USC  Iaage  Processing  Institute  which 
have  been  published  or  accepted  for  publication  during  the  past  six 
aonths,  and  which  have  resulted  fron  ABPA  sponsored  research. 

H.C.  Andrews,  "Numerical  Analysis  Techniques  in  Digital  Inge 
Restoration,"  Proceedings  1975  Syaposiun  on  Circuits  aod  Systems, 
April,  1975. 

H.C.  Andrews,  "HTF  Restoration  by  Pseudoinversion,"  Proceedings  of 
the  International  Optical  Computing  Conference,  April,  1975, 
Washington,  D.C. 

H.C.  Andrews,  Chapter  4,  "Two  Dimensional  Transforms,"  Picture 
Processing  and  Digital  Filtering,  F.S.  Huang,  ed..  Springer  Verlag, 
Hay,  1975. 

H.  C.  Andrews  and  C.  L.  Patterson  "Outer  Product  Expansions  and 
Their  Uses  in  Digital  Image  Processing,"  IEEE  Transactions  on 
Computers,  (accepted  for  publication). 

H.  C.  Andrews  and  C.  L.  Patterson,  "Digital  Interpolation  of 
Discrete  Images,"  IFEE  Transactions  on  Computers  (accepted  for 
publication) . 

H.  C.  Andrews  and  C.  L.  Patterson,  "Singular  Value  Decompositions 
and  Digital  Iaage  Processing,"  IEEE  Transactions  on  Audio,  Speech,  and 


Signal  Processing)  (accepted  for  publication) . 

H.  c.  Andrews  and  B.  R.  Hunt,  Digital  Image  Restoration,  Prentice 
Hall  (accepted  for  publicaton). 

s.  R.  Dashiell  and  A.  A.  Sauchuk,  "Optical  Synthesis  of  Nonlinear 
Nonaonotonic  Functions,"  accepted  for  publication  in  Optics 
Ccaaunicatons. 

w . Frei,  "The  Need  for  a Minimum  Picture  Data  Basis,"  presented  at 
1975  IEEE  Computer  Society  Workshop  on  Machine  Pattern  Analysis,  March 
3-5,  1975. 

W.  Frei,  "Accuracy  Considerations  for  Digitized  laages  and  Hardcopy 
Output,"  presented  at  1975  IEEE  Computer  Society  Workshop  on  Hachine 
Pattern  Analysis,  March  3-5,  1975. 

E.  L.  Hall,  W.  0.  Crawford,  And  F.  E.  Roberts,  "floaent 
Measurements  for  Computer  Texture  Discriaination  in  Chest  X-Rays," 
IEEE  Transactions  Bipaedical  Engineering,  November,  1975. 

E.  L.  Hall,  Z.  H.  Cho,  J.  K.  Chan,  R.  P.  Kruger  and  D.  G. 
McCaughey,  "A  Comparative  Study  of  3-D  Image  Recgnstruction 
Algorithms,"  IEEE  Trans,  on  Nuclear  Science,  March,  1975. 

E.  t.  Hall,  R.  P.  Kruger  and  A.  F.  Turner,  "Automated 
Measurements  frcn  Chest  X-Rays  for  Lanj  Disease  Classification," 
USA-Japan  Computer  Conference,  August,  1975. 


E.  L.  Hall,  R.  P.  Kruger  and  A.  F.  Turner,  "Automated 

* 


Measurements  fron  Chest  X-Rays,"  Proceedings  of  the  Computer 
Applicatons  in  Radiology  Conference,  March,  1975. 

G.  L.  Hall,  S.  B.  Thompson,  G.  Varsi  and  R.  Gaulden,  "Computer 
Measurement  of  Particle  Sizes  in  Electron  Microscope  Images,*1  IEEE 
Trans  on  Systems,  Man  and  Cybernetics,  to  be  published,  1975. 

G.  C.  Huth  and  E.  L.  Hall,  "Coaputer  Tomography  and  its 
Application  to  Nuclear  Medical  Imaging,"  Computers  in  Nuclear 
Medicine,  to  be  published. 

N.  D.  A.  Mascaxenhas  and  9.  K.  Pratt,  "Digital  Image  Restoration 
Under  a Regression  Model,"  IEEE  Transactions  on  Circuits  and  Systems, 
March,  1975. 

N.  E.  Nahi  and  N.  Naraghi,  "A  General  Image  Estimation  Algorithm 
Applicable  to  Multiplicative  and  Non-Gaussian  Noise,"  Proceedings  of 
18th  Midvest  Symposium  on  Circuits  and  Systems,  August  11-12,  1975, 
Concorshia  Univ.,  Montreal  P.Q.,  Canada. 

N.  E.  Nahi  and  A.  Habibi,  "Nonlinear  Recursive  Image  Enhancement," 
IEEE  Transactions  on  Circuits  and  Systems,  March,  1975. 

R.  Nevatia,  T.  0.  Binford,  "Recognition  and  Description  of  Cpuplex 


Curved  Objects," 

Fifth 

Annual 

Symposium 

on  Imagery  Pattern 

Recognition,  University  of 

Maryland, 

April  17-18, 

1975. 

R.  Navatia  and  1. 

0. 

Binford, 

"Recognition 

and  Description  of 

Complex  Curved  Objects",  Fifth  Annual  Symposium  on  Automatic  Imagery 
Pattern  Recognition,  Univ.  of  Maryland,  April  1975. 
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B.  Nevatia,  "Object  Boundary  Deteraination  in  a Tejtured 
Environment,"  (Tc  be  presented)  Annual  ACM  Conference,  Minneapolis, 

October  1975. 

B.  Nevatia,  "Eepth  Heasureaent  by  Motion  Stereo,"  Accepted  for 
publication  in  Ccaputer  Graphics  and  Iaage  Processing. 

B.  Mevatia,  "Structured  Descriptions  cf  Coupler  Curved  Objects  for 
Recognition  and  Visual  Heaory,"  Accepted  for  publication  as  a bpok  by 
Birkhauser-Verlag,  Basle,  Switzerland. 

M . R.  Pratt  and  M.  Uuhns,  "DPCN  Quantization  Error  Reduction  for 
Iuage  Coding,"  Society  of  Photo-Optical  Instruaentation  Engineers, 

19th  Annual  Technical  Symposium,  San  Diego,  August,  1975. 

H.  K.  Pratt  and  C.  E.  Hancill,  "Spectral  Estivation  Technigues  for 
the  Spectral  Calibration  of  a Color  Iaage  Scanner,"  Applied  Cptics,  j 

Noveaber,  1975. 

M.  K.  Pratt,  "Vector  Space  Foruulation  of  Two  Diuensional  Signal 
Processing  Operations,  Journal  Coaputer  Graphics  and  Iaage  processing, 

Acadenic  Press,  March,  1S75. 

J.  A.  Roese,  W.  K.  Pratt,  G.  S.  Robinson,  A.  Habibi, 

"Interfraae  Transfora  Coding  and  Predictive  Coding  Metods,"  IEEE 
Internatonal  Couaunicatons  Conference,  San  Francisco,  June,  1975. 

J.  A.  Roese,  G.  S.  Robinson,  "Coabined  Spatial  and  Temporal  Coding 
of  Iaage  Seguences",  19th  Annual  SPEI  Syaposiun  on  Efficient 
Transmission  of  Pictorial  Information,  San  Diego,  Calif.,  August 
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21-22,  1975. 


A.  A.  Sawchuk  and  M.  J.  Peyrovian,  "Restoration  of  Astigmatism  and 
Curvature  of  Field",  Journal  of  the  Optical  Society  of  America,  vol. 
65,  1975. 

A.  A.  Sawchuk  and  S.  R.  Oashiell,  "Nonmonotonic  Nonlinearities  in 
Optical  Processing",  Proc.  IEES  International  Optical  Computing 
Conference,  Washington,  D.C.,  April  23-25,  1975. 

N.  s.  Thompson,  A.  F.  Turner,  and  R.  P.  Kruger,  "Automated  Chest 
Radographic  Diagnosis,"  accepted  for  publication.  Investigative 
Radiology. 

V.  E.  Thompson,  A.  I.  Zobrist,  "Building  a Distance  Punction  for 
Gestalt  Grouping,"  accepted  for  publication  IEEE  Transactions  on 
Computers. 


3.  Image  Coding  Projects 


The  effort  in  iaage  coding  is  directed  toward  the  research  and 
development  of  iaage  coding  systems  providing  a transmission  bit  rate 
reduction  and  tolerance  to  channel  errors.  Coding  systems  are  under 
investigation  for:  monochrome  and  color  imagery;  slow  scan  and  real 
tine  television;  and  information  preserving  and  controlled  fidelity 
operation.  Results  of  this  research  study  during  the  past  six  months 
are  summarized  here  and  presented  in  detail  in  subsequent  sections. 

3.1  Singular  Value  Decomposition  Image  Coding 
Harry  C.  Andrews 

The  singular  value  decomposition  algorithm  (also  referred  to  as 
SVD)  is  a computational  algorithm  developed  for  a variety  of 
applications  including  matrix  diagonalization,  regression,  and 
pseudoinversion  £1,2].  The  algorithm  has  also  been  suggested  for 
image  coding  £3,4].  By  approaching  the  image  coding  task  from  a 
viewpoint  of  numerical  analysis,  it  is  possible  to  formulate  a 
solution  in  terms  of  least  squares  methods  which  results  in 
deterministically  best  truncation  errors  over  all  other  unitary 
transforms  [6].  A discrete  iaage  may  be  considered  to  be  an  array  of 
nonnegative  scalar  entries  formed  into  a matrix.  Let  such  an  iaage 
matrix  be  designated  as  G.  Without  less  of  generality,  let  G be 
square  with  a quadratic  fora  given  by 


(1) 


G 


AaB 


T 


where  the  A_  and  B_  matrices  are  assumed  unitary.  Salving  for  a it  is 
observe)  that 


a = 


T 

AG  B 


(2) 


The  £ matrix  is  seen  to  be  the  "transform**  of  the  iaage  matrix  where 
i transforms  the  coluans  of  the  iaage  and  JL  transforas  the  rows  of 
the  iaage.  A list  of  traditional  transform  techniques  is  presented  in 
Table  1 indicating  sane  of  the  properties  of  the  individual  transfora 
methods.  The  entries  are  listed  in  terns  of  general  decreasing 
usefulness  as  decorrelation  devices  as  well  as  decreasing  complexity. 

The  first  entry  in  the  table  is  the  one  of  interest  here  and  has 
decidedly  different  transform  properties  froa  the  remaining.  The 
singular  value  decomposition  (SVD)  method  has  the  unique  property  that 
the  coefficient  aatrix  A_  , is  diagonal  with  at  most  only  N nonzero 
entries.  The  definition  of  this  transform  is  given  by 

a=  A*  = UTG  V (3) 

where 
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and 


G 


V A V 


(**b) 


The  A matrix  is  diagonal  and  comprises  the  singular  values  of  the 

picture  matrix  G,  while  JJ,  and  ace  the  respective  singular  vector 
T T 

matrices  of  GG  and  G G#  and  are  orthogonal  as  a result  of  the 

T T 

nonnegative  definiteness  cf  GG_  and  G G.  Because  of  these  properties 
of  0 and  V it  is  possible  to  solve  for  the  image  matrix  such  that 


k T 

G = U A 2 V 


(5a) 


or  equivalently 


G = V"'  A.2  u.vT  (5b) 

— x — i~i 

i=  1 

Where  H ^ N and  represents  the  rank  or  number  of  nonzero  singular 
values  A.  . The  coding  implications  are  that  one  must  trausmit  the  2N 
singular  vectors  as  veil  as  N singular  values  for  image  reconstruction 
at  the  receiver.  Figure  1 presents  a pictorial  representation  of  the 
singular  value  decomposition.  Traditional  image  transform  methods 
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usually  break  an  image  up  into  subblocks  for  ease  in  hardware 
implementations.  This  technique  is  developed  here  because  the 
computational  expenses  for  large  image  singular  value  deepapositions 

is  great.  Specifically,  if  an  N x K image  is  broken  into  N x N 

3 

sutblocks,  then  each  subblock  takes  on  the  order  of  N computations  to 

z 

get  to  the  SVD  domain.  Since  there  are  (M/N)  such  subblocks,  a total 

7 3 

of  Mm  computations  are  required  as  compared  to  H computations  if  the 

entire  image  mere  decomposed.  A similar  comparison  exists  for  fast 

2 

computational  transforms  which  require  K log  H total  subblock 

operations  for  the  « x M image.  Thus  the  number  of  computations  for 

2 2 , 

SVD  compared  to  fast  transforms  is  M N vs  M log  N.  The  ratio  of  N/log 
N increase  to  implement  the  SVD  transform  on  16  x 16  subblock  sizes  is 
only  a factor  of  4 for  SVD  versus  Fourier,  cosine,  Walsh,  or  the  like. 

Figure  2 contains  a block  diagram  of  the  SVD  ceding  scheme.  The 
major  components  at  the  transmitter  consist  of  the  SVD  domain 

transformer,  a possible  truncator,  and  parallel  singular  value  and 
singular  vector  coders.  The  SVD  transformer,  as  discussed  above, 
would  require  on  the  order  of  four  times  the  number  of  real 

computations  compared  to  a real  NZlog  N transform  algorithm.  The 

truncator  is  included  in  the  diagram  to  emphasize  the  tremendously 
large  energy  compaction  property  of  the  SVD  technique.  From  eg.  (6) 
the  truncated  image  G may  contain  an  extremely  large  amount  of 

ociginal  image  energy  in  a very  few  nuxber  of  singular  values. 

The  two  remaining  blocks  in  the  coder  concern  themselves  with  the 
singular  value  coding  and  singular  vector  coding.  In  the  former  the 


Figure  3. 1-2.  S S/D  Coding  Scheme 


large  dynamic  range  of  the  singular  values  indicates  a certain  amount 

of  care  needed  for  coding,  but  because  there  are  so  feu  large  dynamic 

2 

ranged  coefficients,  (N  vs.N  ) the  total  bit  requirement  still  remains 
low.  The  singular  vector  coding  algorithm  is  broken  into  two 
components,  that  associated  with  the  first  singular  plane  (or 
eigenimage)  and  that  associated  with  the  remaining  eigenimage  planes. 
Because  each  of  these  planes  (actually  two  vectors  which  when,  outer 
producted,  produce  a plane)  is  orthonormal,  the  scalar  entries  in  the 
singular  vectors  are  quite  well  behaved,  and  lend  themselves  to  easy 
requantization. 

Using  the  basic  configuration  of  figure  2,  the  number  of  bits 
necessary  for  transmission  of  a subblock  then  becomes  a function  of 
the  truncation,  if  any,  the  bit  distribution  over  the  singular  values, 
and  the  bit  distribution  over  the  singular  vectors.  Typical 
distributions  on  the  singular  values  track  the  variance  of  these 
values,  and,  in  fact,  tend  to  be  proportional  to  the  value  itself. 
For  the  singular  vectors,  two  more  bits  are  provided  for  scalar 
entries  of  the  first  eigenimage  than  for  subsequent  eigenimages.  in 
addition,  because  the  singular  vectors  are  orthonormal,  one  need  not 
transmit  N scalar  values  per  vector  but  only  N-i-1  such  valves  for  the 
i-th  vector.  (Orthcnocmality  reduces  the  degrees  of  freedom  cn  the 
vectors  such  that  the  vectors  can  be  completed  at  the  receiver.) 

In  order  to  develop  efficient  quantizers  and  coders  for  the  SVD 
domain,  a test  image  (512  x 512)  was  broken  into  16  x 16  subblocks  and 
data  was  gathered  over  each  subblock  of  the  entire  frame.  Statistics 
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describing  the  singular  values  and  vectors  were  then  gathered  and  are 
described  here.  For  16  x 16  subblocks,  one  obtains  16  singular  values 
and  raonotonic  decreasingly  ordered  means  and  variances  of  each 
singular  value  can  be  computed.  The  exceedingly  large  dynamic  range 
of  between  4 and  5 orders  cf  magnitude  indicates  the  need  fpr  variable 
bit  coding  as  a function  of  the  singular  value  index.  The 
distribution  of  the  singular  values  naturally  is  one  sided  (no 
negative  entries)  and  appears  as  a curve  intermediate  between  a broad 
Gaussian  and  uniform  density. 

The  statistics  describing  the  singular  vectors  are  mmch  better 
behaved.  Figure  3 presents  two  specific  singular  vectors  from  a 
particular  suttlock  as  an  illustration  of  the  shape  of  these 
parameters.  The  singular  vectors  tend  to  be  well  behaved  in  their 
range  and  tend  tc  have  an  increasing  number  of  zero  crossings  as  a 
function  of  increasing  index.  In  fact  it  is  known  that  the  first 
singular  vector  never  exhibits  zero  crossings  when  the  smbblock  is 
ncnnegative  (as  it  always  is  for  imagery)  [7].  Thus  the  lower  indexed 
vectors  tend  to  have  a great  deal  of  adjacent  sample  correlation. 

Since  the  first  vector  for  both  (Land  J l are  guaranteed  to  have  no 

zero  crossings  (similar  to  the  dc  vectors  of  Fourier,  Walsh,  cosine, 

etc.  transforms),  these  vectors  form  a separate  set  of  statistical 

parameters  from  the  remaining  set.  The  mean  vector  over  all  subblocks 

in  the  test  image  bedomes  a constant  value  of  0.25  with  a very  tight 

-3 

variation  provided  by  a variance  of  10  . Naturally  a given  first 
singular  vector  will  not,  in  fact,  be  a constant  of  value  0.25  (the 
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appropriate  notaalized  value  to  guarantee  or thonoraality) , but  Hill 

have  variations  which  when  weighted  by  the  corresponding  large 

singular  value  will  appear  guite  different  from  a constant.  In  fact 

the  distribution  pf  the  scalar  values  defining  the  entries  in  the 

first  singular  vector  are  very  close  to  Gaussian  with  parameters 
-3 

N (0.  25,10  ). 

The  renaining  eigenvectors  are  also  guite  well  behaved  with  the 
average  or  aean  of  each  cf  these  vectors  being  the  zero  vector.  The 
variance  of  these  singular  vectors  is  on  the  order  of  10  1 and  the 
scalar  entries  which  comprise  these  vectors  are  also  close  to  ncrnally 
distributed  with  *((0,10  *)•  Because  of  the  difference  in  the 

statistics  of  the  first  singular  vector  with  those  of  the  renaining 
singular  vectors,  they  are  coded  separately  as  indicated  in  the  block 
diagran  of  figure  2. 

one  inage  is  used  fcr  experimental  purposes  here.  its  SVO 
structure  is  revealed  in  figure  4.  In  figure  4 the  inage  is  broken 
into  16  x 16  sutblocks  and  each  subblock  is  deconposed  into  its  16 
singular  values  and  associated  32  singular  vectors.  The  first, 
second,  third,  and  fourth  eigeninages  obtained  frcn  the  corresponding 
singular  vectors  are  displayed  in  the  figure.  The  first  plane  has  no 
zero  crossings  and  conseguently  the  display  of  negative  data  is  not 
necessary.  In  the  three  retaining  planes,  a linear  dislay  is 
presented  with  negative  values  being  dark  and  positive  values  being 
light.  Notice  that  considerable  recognition  of  the  original  scene  is 
available  in  the  first  eigeninage  and  subseguent  eigeninages  tend  to 
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provide  more  and  more  zeto  crossing  information  which  can  he  likened 
to  higher  frequency  information. 

k complete  STD  coding  system  has  been  simulated  according  to  the 
block  diagraa  cf  figure  2.  The  ceding  strategy  used  a linear  PCM 
quantizer  with  variable  bit  assignment  on  the  singular  valmes  and  a 
Nax  quantizer  £5]  with  variable  bit  assignment  on  the  PCM  values  of 
the  entries  in  the  singular  vectors.  A variety  of  bit  assignments 
were  investigated,  and  an  optimization  routine  in  terms  of  mean  square 
error  measured  the  best  kit  assignment.  Figure  5 provides  some 
performance  curves  developed  during  the  optimization  process.  The  two 
lower  curves  indicate  the  truncation  effects  as  the  nuirber  of 
equivalent  bits  per  pixel  are  increased.  The  uppermost  curve 
illustrates  the  lean  square  error  using  a linear  quantizer  on  the 
singular  values.  The  Max  quantizer  curve  indicates  about  a 0.20%  mean 
square  error  improvement  over  the  linear  curve  and  is  only  about  d.20% 
worse  (or  introduces  0.20*  more  mean  square  error)  than  the  truncated 
but  uncoded  curves.  Pictorial  results,  from  which  the  upper  two 
curves  are  derived,  are  presented  in  figure  6.  Here  the  percentage 
mean  square  errors  and  bit  rates  per  pixsl  are  listed  under  the 
respective  coded  images  for  both  linear  and  Max  quantization  on  the 
singular  vectors. 

In  concluding  this  section  it  is  important  to  emphasize  a few 
points.  First,  the  work  is  incomplete  and  it  is  premature  to  base 
any  conclusions  cn  the  viability  of  STD  coding  in  competition  with 
ether  existing  technigues.  It  is  fair  to  say  that  if  as  much  effort 
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BITS  PER  PIXEL  (bpp) 

e mean  square  error  versus  hits  per  pixel. 


Figure  3.1-6.  "Baboon"  Image  with  SVD  coding  using  "Baboon"  statistics 


is  put  into  investi gatiu g the  potential  for  SVD  coding  as  has  been  put 
into  traditional  transfora  methods,  then  considerable  iaproveaent  over 


the  results  presented  here  can  be  expected.  However,  algorithmic 
implementation  sight  become  quite  cosplex.  On  the  other  hand  five 
years  ago  realtime  (video  bandwidth  rate)  FPT  transfora  cpders  were 
thought  to  be  too  complex,  and  yet  they  exist  today.  Consequently 
only  tiae  and  future  study  will  tell  whether  SVD  coding  becoaes  a 
practical  reality. 
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3.2  Restoration  for  Binary  Syaaetric  Channel  Errors 
Michael  N . Huhns 

A previous  report  [1]  has  presepted  and  analyzed  a technique  for 
restoring  the  output  of  a quantizer  so  that  the  result  more  accurately 
Batches  the  quantizers  input  with  respect  to  a mean-square  error 
criterion.  The  restoration  is  obtained  by  the  use  of 


E {x  | x e R] 


/ x p(x)dx 
^ R 

/ p(x)dx 


(1)  I 

i 


whare  R is  a region  in  N-space  to  which  an  N x 1 vector  x is  assigned 
during  quantization,  and  p(x)  is  the  multidimensional  probability 
density  function  of  x.  The  restoration  is  based  essentially  upon 
exact  knowledge  of  the  quantizer  output.  A similar,  but  more 
difficult  problem  results  when  the  quantizer  output  is  not  known 
exactly.  This  could  occur,  for  example,  when  the  quantizer  output  is 
transmitted  over  a npisy  channel.  The  first  section  in  this  report 
explores  the  effect  of  channel  errors  on  the  restorations  obtained 
using  eq.(1).  The  next  section  examines  a technique  that 
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statistically  compensates  for  the  effect  of  channel  errors. 

Effects  of  Channel  Errors  on  Quantized  Signals:  In  this  analysis, 
channel  errors  are  assumed  to  be  modelled  by  a binary  symmetric 
channel  (BSC)  [2].  The  characteristics  of  this  type  of  channel  are 
shown  in  figure  1.  The  channel  is  discrete  and  memoryless  and  can  be 
specified  by  a transition  probability  assignment  F(j|k),  for  j,k=0,1, 
as 


Since  the  channel  is  memoryless,  the  probability  of  an  output  sequence 
2=  (Zj  ,z2,  . . . ,zN)  , given  an  input  sequence  jt=  (Xj  , *2,  ...  .x^) , is  given 
by 

N 

p(z  | x)  =TT  p(z.  |x.)  (3) 

i=  1 11 


, Based  on  this  definition,  a BSC  was  computer  simulated  with  the 

channel  error  probability,  p,  chosen  to  be  0.01.  The  simulated 
channel  was  then  applied  tc  transform  coded  images.  Three  images  were 
zonal  transfer*  coded  in  16  x If  blocks  and  their  quantized  transform 

‘ domain  components  were  encoded  by  assigning  each  a binary  code  word. 

* . 

« The  resulting  sequence  cf  binary  digits  was  operated  on  by  the 

* 


* 

* 


1 


simulated  channel.  The  er lor-corrupted  bit  streaa  was  then  either 
decoded  directly,  as  shown  in  figures  2a,  2c,  and  2e,  or  restored  by 
the  use  of  eq.(1)  to  reduce  the  effects  of  the  quantization  process. 
Figure  3 contains  a scheaatic  of  this  procedure.  The  decoded  iaages 
with  the  quantization  effects  reduced  are  shown  in  figures  2b,  2d,  and 
2f. 

Bit  errors  in  transfora  coding  that  arise  due  to  a binary 
syaaetric  channel  are  seen  to  result  in  an  emphasis  of  the  block 
structure  and  a subjective  error  that  extends  over  the  entire  block. 
This  latter  effect  occurs  because  inverse  transforaing  a block 
containing  an  error  distributes  this  error  over  all  the  resultant 
iaage  doaain  components.  The  reconstruction  technique  iaplied  by  eq. 
(1)  is  thus  insensitive  to  channel  errors.  Since  it  provides  visual 
and  aean-square  erxor  inproveaents  in  noise-free  cases,  it  can  be 
utilized  equally  sell  in  noisy  enviconaents. 

Reconstruction  of  Quantized  and  Transaitted  Signals:  The  previous 
section  deaonstrated  that  channel  errors  do  not  adversely  affect  the 
performance  of  the  restoration  technique  derived  previously.  However, 
this  technique  does  nothing  to  aaeliorate  the  effects  of  the  channel 
errors.  This  is  because  the  fundamental  restoration  foraula  presented 
in  eg.  (1)  was  derived  without  any  consideration  of  channel  structure. 
By  including  the  channel  structure  in  the  derivation,  the  resultant 
restoration  technique  can  simultaneously  reduce  the  effects  of  the 
quantization  process  and  aitigate  the  effects  of  channel  errors. 
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(c)  Quantized  0.  5 bit/pixel  (d)  Restored  0.  5 bit/pixel 

P =0. 01  P =0.  01 

e e 


(e)  Quantized  0.  5 bit/pixel  (f)  Restored  0.  5 bit/pixel 

P =0.01  * P =0.01 

e e 


Figure  3.  2-2,  Minimum  mean  square  error  restoration  of  Hadamard 
transformed  zonal  quantized  images. 


-28-  j 


DATA 

SOURCE 


QUANTIZER 


DATA 

ENCODER 


Let  the  output  cf  a data  source  (this  output  could  consist  of 
DPCN  samples,  PCM  samples,  or  transform  domain  samples)  he  denoted  by 
— = (*  i#x  2*  ’ * * *xn)  an<1  <JescE  a probability  density  function  p (x) . 
The  ceconstructicn  of  x_,  after  x has  been  quantized  to  one  of  M 
regions  and  channel-error  corrupted,  is  denoted  by  z= (z^ ,22»* *•  ,2^) 
for  k=1,2,...,M  (refer  to  figure  3).  The  mean-square  error  that 
results  from  this  process  is 


(4) 


This  error  can  be  minimized  by  proper  choice  of  the  restoration 
points,  z^.  Setting  the  partial  derivatives  of  this  error  vith 
respect  to  zk  equal  to  zero  yields 


M 


m=  1 

Zk  = M 


Y'  p(m|k)  f xp(x)dx 


(5) 


yi  p(m|k)^  p(x)dx 
m=  1 


for  k=1,2,...,M.  This  expression  is  the  noisy  channel  version  of 
eq.(1)  and  provides  a minimum  mean-square  error  estimate  of  the  input 
to  a quantizer  based  on  the  output  of  a noisy  channel,  the 
characteristics  of  the  quantizer,  and  the  a priori  statistics  pf  the 
input.  This  equation  is  also  a multidimensional  version  of  a result 
first  derived  in  [3].  For  a noiseless  channel,  the  channel  matrix  P 
becomes  the  identity  matrix  and  eq.  (5)  reduces  to  eq. (1) . When  the 
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pcobability  vclaae  integrals  in  the  denominator  of  eg.  (5)  ace  all 
equal,  which  is  approximately  true  foe  Max  quantization,  the 
restocation  equation  simplifies  to 


I xp(x)dx 
M *11 

zk  = S 

m=  1 / p(x)dx_ 

•'n 


= *2  p(mlk)yt 


where  ymis  given  by  eq.  (1).  This  result  holds  for  maximum  output 
entropy  quantizers  and  two-level  symmetrical  quantizers,  and  is 
approximately  correct  for  many  other  types. 

A signal  that  has  been  quantized  and  then  transmitted  over  a 
noisy  channel  can  thus  be  optimally  restored  by  utilizing  eg. (5).  The 
restoration  scluticns  found  earlier  for  Gaussian  and  Laplacian 
probability  density  functions  (see  [4]  an  1 [5],  respectively)  can  be 
substituted  directly  into  eg.  (5)  once  the  transition  matrix  for  the 
channel  has  been  determined.  The  resultant  estimator  can  then  be  used 
to  restore  the  outputs  of  transform  and  DPCJ1  coders  that  have  been 
degrade!  by  channel  errors. 
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3.3  Interframe  Iaage  Coding 

Guner  S.  Robinson  and  John  A.  Roese 

Interfcame  coding  of  digital  iaage  sequences  encompasses  those 
techniques  which  make  use  ◦£  the  high  correlation  that  exists  between 
pixel  amplitudes  in  successive  frames.  Intraframe  coding  techniques 
that  exploit  spatial  correlations  can,  in  principle,  be  extended  to 
♦ This  research  is  partially  supported  by  the  Naval  Undersea  Center, 
San  Ciejo,  California. 


include  correlations  in  the  temporal  dona  in.  Previous  research  in  the 
area  of  three-diaeosional  Pourier  and  Hadaaard  transformations  has 
indicated  that  bit  rates  can  be  reduced  by  a factor  of  five  by 
incorporating  correlations  in  the  teaporal  direction  [1].  Hoaever, 
three-dimensional  transfoca  systems  are  unattractive  since  they 
require  large  amounts  of  data  storage  and  excessive  computation. 

To  alleviate  the  problems  associated  vith  three-dimensional 
transform  systems,  nev  hybrid  (two-diaensional  transform) /DPCH  iaage 
coding  systems  have  been  developed  [2].  These  systems  utilize  both 
spatial  and  teaporal  correlations  while  greatly  reducing  aeaory 
storage  and  coaputat ional  requirements.  A block  diagram  for  a hybrid 
(two-dimensional  transf or  a) /QPCN  system  is  shown  as  figure  1.  In 
present  implementations  of  this  system,  either  a two-diaensional 
cosine  or  Fourier  transformation  is  performed  on  16  x 16  subblocks. 
DFCM  linear  predictive  coding  in  the  temporal  domain  is  then  applied 
to  the  transform  coefficients  of  each  subblock.  For  notational 
convenience,  the  hybrid  interframe  coders  employing  two-dimensional 
Fourier  transforms  will  be  denoted  as  FFD  and  those  using 
two-d imensicnal  cosine  transforms  as  CCD.  The  FFD  and  CCD  coders  are 


adaptive  in  the  sense  that  statistics  of  the  transform  coefficient 
differences  of  each  subblcck  are  computed  prior  to  encoding  the 
transform  coefficients  in  the  temporal  direction  by  parallel  DPCH 
coders.  At  the  receiver,  the  transmitted  transform  coefficients  are 
decoded  and  a replica  of  each  frame  is  reconstructed  by  the 
appropriate  inverse  two-dimensional  transformation.  These  systems 
require  only  a single  frame  of  storage  and  involve  significantly  less 


f 


? 


memory  and  fewer  computations  than  three-diaensional  transfora  coding 
technigues. 

Operational  Bodes:  At  least  three  operational  nodes  have  been 
identified  for  the  hybrid  interfrane  coding  systeas.  These 
operational  nodes  depend  on  the  initial  conditions  assuaed  for  the 
previous  coder.  The  initial  conditions  are: 

a.  No  apriori  inforaation  available  at  the  receiver; 

b.  Liaited  infpraaticn  (such  as  mean,  variance  and  teaporal 

correlations  based  on  a statistical  aodel) 
available  at  the  recaiver;  and 

c.  First  fraae  available  at  the  receiver. 

In  the  no  apriori  inforaation  available  case,  several  frames  are 
reguired  for  the  hybrid  coder  to  settle.  However,  it  has  been 
experimentally  verified  that  in  the  reaaining  two  cases,  nearly  stable 
coder  performance  is  achieved  within  the  first  U to  6 franes.  Proa 
operational  considerations,  the  third  set  of  initial  conditions  is  the 
most  realistic  as  periodic  full  frane  updating  will  be  reguired  to 
eliminate  the  cuvulative  effects  of  channel  noise. 

Mathematical  Formulation:  Let  f (x,y)  denote  a two-dimensional 
array  of  intensity  values  on  an  N x N subblock  of  a digital  television 
image  of  size  M x M.  Typical  values  for  H and  N are  256  and  16, 
respectively.  Let  F(u,v)  be  the  two-dimensional  array  obtained  by 
taking  the  two-dimensional  transform  of  f (x,y) . In  the  case  of  the 
two-dimensional  discrete  Fourier  transform,  the  expressions  relating 
f (x, y)  and  F (u, v)  are 


-34- 


(1) 


N-l  N-l 

F(u,v)  =~  ^2  ^ f(x,  y)  exp  (ux  + vy)l 

N x=0  y- 0 L J 


and 


N-l  N-l 

f(x,  y)  E v)exp  1+  (ux  + vy)l  (2) 

u=0v--0 

foe  u, v, x, y«0,  1 , . . . ,h-  1.  Por  image  processing  applications,  f (x,y)  is 

a positive  real  function  representing  brightness  of  the  spatial 

sample.  The  two-dimensional  Fourier  tnnsform  of  a real-valued 

function  has  the  conjugate  symmetry  property.  Also,  the  Fpurier 

2 

transfora  consists  of  2M  components,  i.e.,  the  real  and  imaginary  or 

magnitude  and  phase  components  of  each  spatial  frequency.  However,  as 

2 

a result  of  the  conjugate  eymmetry  properties  mentioned  above,  only  N 
components  are  reguired  to  completely  define  the  Fourier  transform 
[3]. 


In  the  case  cl  the  Fourier  transform,  a shift  in  the 
spatial-domain  variables  results  in  a multiplication  of  the  Fourier 
transform  of  the  on-shifted  image  by  a phase  factor.  If  the  input 
imge  f(x,y,t1)  is  shifted  by  the  amount  xQ  in  the  x-direction  and  yQ 
in  the  y-direction  between  times  t ^ and  t2,  then  the  Fourier  transform 
of  the  shifted  image  is  given  by 

F(u,v,t2)  = Ffu.v.tj)  exp  (uxQ  + vyQ)j  (3) 
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This  shifting  property  is  expected  to  be  useful  in  detecting  and 
compensating  for  effects  cf  motion  between  frames  since  many  types  of 
motion,  such  as  panned  moticn,  produce  significant  changes  in  phase 
components  but  small  changes  in  amplitude  components.  Thus, 
compensation  for  camera  platform  moticn  could  be  implemented  directly 
in  the  array  of  phase  eexponents  by  application  of  appropriate  phase 
correction  factors. 

The  two-dimensicnal  Fourier  transform  ?{u,v)  of  a spatial  signal 
function  f(x,y)  is  separable,  i.e.,  it  can  be  computed  as  two 
sequential  one-dimensional  transforms  since  the  Pourier  kernel,  is 
separable  and  symmetric.  Thus,  the  basic  one-diaensional  discrete 
Fourier  kernel  transform  that  must  be  performed  is 


*■<»>  =|rX  ,(x)  e*p  (-t?  ”*) 


for  us0f 

in  the  case  of  the  discrete  Cosine  transform,  the  one-dimensional 


transform  is 


F(u)  2 f(x)  cos  ((2VN1)U") 

x=  0 v ' 


for  u=0, 1, . . . ,N- 1.  The  cosine  transform  is  also  separable  and  a 

two-dimensional  discrete  cosine  transform  of  an  N x H subblcck  results 

2 . . 

in  N real  coefficients. 
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Experimental  evidence  derived  from  transmission  of  a typical 
"head  and  shoulders"  picture  telephone  scene  has  shown  that  the  frame 
difference  signal  has  a probability  density  closely  approximated  by  a 


r 


double  sided  exponential  function  [4].  The  optimum  minimum  mean 
square  error  quantizer  for  this  distribution  has  been  found  tc  be  a 
uniform  quantizer  combined  with  a companding  of  the  frame  difference 
signal  £ 5 ]. 


Since  the  variances  of  the  transform  domain  coefficient 
differences  are  different,  it  is  necessary  to  use  different  quantizer 
parameters  foe  each  coefficient..  Each  coefficient  difference  signal 
is  allocated  a number  of  bits  proportional  to  the  estimated  variance 
in  accordance  with  an  optimum  bit  assignment  algorithm. 


Fidelity  Criteria:  In  figure  1,  differences  between  the  input 

A 

signal  f(x,y,t)  and  output  signal  f(x,y,t)  are  due  to  two  sources: 
quantization  errors  and  channel  noise  errors.  To  evaluate  Coding 
efficiency  of  the  hybrid  encoders,  two  objective  criteria  were  used. 
The  first  criterion,  NHSE,  is  a measure  of  the  mean  sgaare  error 
between  f(x,y,t)  and  f(x,y,t)  averaged  over  an  entire  frame  of  size  N 
x r.  Normalization  is  achieved  by  dividing  the  mean  square  error  by 
the  mean  signal  energy  within  the  frame  to  give 


M-l  M-l 
2 

NMSE  = — — ~ ~~Y 
M 


1V1-1  AVI-1 

-4  EE  ff(x,y,t)  - f(x,y,t)l 

M x=0  v=0  L J 

1 M-l  M-l  r i 2 

/ EE  [«*■?•*>] 


x=0  x=0 


(6) 
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HYBRID  (TWO-DIMENSIONAL  TRANSFORM)/DPCM  CODER 


1 


The  second  criterion,  SNR , measures  the  ratio  of  peak-to-peak  signal 
to  RMS  noise  as  defined  by 


SNR  = -10  log1Q 


M-1M-1 

EE  [ f(x,y,t) 

x=  0 x=  0 


255 


(7) 
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Figures  2a  and  2b  are  graphs  illustrating  the  coding  efficiency 
of  the  hybrid  FFD  and  CCD  coders  at  various  bit  rates  in  the  interval 
0.1  to  1.0  bits/pixel/fraae.  To  perform  this  series  of  experiments,  a 
256  x 256  resolution  data  base  consisting  of  16  consecutive  fraaes  of 
a 24  fraaes  per  second  (fps)  aotion  picture  was  digitized.  /Initial 
conditions  assuaed  were  that  the  first  frame  was  available  at  the 
receiver. 


Photographs  pf  fraae  nuaber  16  after  coding  by  the  FFD  and  CCD 
coders  at  average  final  bit  rates  of  1.0,  0.5,  0.25,  and  0.1  are  shown 
as  figures  3 and  4.  The  results  shown  in  figure  3 for  the  FFD  coder 
were  obtained  by  Coding  the  real  and  imaginary  components  pf  the 
Fourier  coefficients  by  assigning  half  of  the  available  bits  to  each 
coaponent. 


Noise  Immunity;  Performance  of  the  FFD  and  CCD  hybrid  interfraae 
coders  was  investigated  in  the  presence  of  channel  noise.  In  order  to 
study  the  effect  of  channel  noise,  a binary  symmetric  channel  was 
siaulated.  The  channel  is  assuaed  to  operate  on  each  binary  digit 
independently,  changing  each  digit  frca  0 to  1 or  froa  1 to  0 with 


5£aaMBBS2-^5 
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(a)  1.  0 bit  a/pixel/frame 


(b)  0.  5 bits /pixel/frame 


Figure  3.  3-3  . FFD  coder  for  frame  16. 
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(a)  1.  0 bits /pixel /frame 


(b)  0.  5 bits ^>ixel /frame 
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probability  p and  leaving  the  digit  unchanged  with  probability  1-p. 
At  the  receiver,  the  encoded  picture  is  reconstructed  from  the  string 
of  binary  digits,  including  errors,  transmitted  across  the  channel. 


4 


Degradations  due  to  channel  noise  probabilities,  p of  zero,  10 
-2 

and  10  for  the  FFD  and  CCD  coders  at  average  bit  rates  of  1.0  and 
0.25  bits/pixel/fraae  are  shown  in  figures  5 and  6.  The  generally 
monotonically  increasing  character  cf  these  curves  illustrates  the 
fact  that  once  an  error  has  occurred,  it  tends  to  propagate  in  the 
temporal  direction  until  corrected  by  a frame  refresh. 

Resulting  pictures  shew  that,  fer  both  coder  implementations 
studied,  minimal  image  degradation  occurred  for  channel  error 

_3 

probability  cf  10  cr  less. 

Photographs  corresponding  to  average  bit  rates  of  1.0  and  0.25 

bits/pixel/frame  for  the  FFD  and  CCD  coders  with  channel  error 

-3  -2 

probabilities  of  10  and  1C  are  shown  in  figures  7 and  3. 

Bit  Transfer  Rate;  In  keeping  with  the  previously  mentioned 
objective  of  ninimizing  the  number  of  bits  transmitted  while  retaining 
image  fidelity,  a series  of  experiments  was  performed  in  which  certain 
bit  transfer  rates  (BTR)  across  the  channel  were  fixed.  The  ETR  is 
defined  as  the  pteduct  of  average  pixel  bit  rate  per  frame  and  frame 
rate  and  has  units  of  bits/pixal/sec. 

The  available  16  frame  test  data  base  was  extracted  from  a 24  tps 
motion  picture  sejuence.  dy  employing  frame  skipping  techniques. 
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(a)  1.0BITS/PIXEUFRAME 


(b)  0 25  BITS/PIXEL/FRAME 
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Figure  3.3-5.  Effects  of  Channel  Noise  For  Fourier/Fourier ^PCM  Coder 
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(c)  0.  25  bits^rixel/frame 

_ 3 

p = 10 


(d)  0.25  bit s /pixel /frame 

p=  io"2 


Figure  3.  3-7.  FFD  coder  with  channel  noise. 
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(a)  1.  0 bits  /pixel /frame 
p=10"3 


(b)  1.0  bits  /pixel/frame 

p = io~2 


p Figure  3.  3-8.  CCD  coder  with  channel  noise. 


temporal  subsampling  was  used  to  simulate  short  12,  8 and  6 fps 
sequences  from  the  16  frame  test  data  base. 

Average  bit  rates  in  the  interval  0.083  to  1.333  bits/pixel/fraae 
were  used  in  ccajunction  kith  the  four  fraae  rates  mentioned  above  to 
perfcra  siaulaticns  with  BT R values  of  8,  6,  4 and  2 bits/pixel/sec. 
Results  of  these  experiments  for  4 bits/pixel/sec  are  shown  in  figure 
9.  For  all  cases  exaained,  the  graphs  show  that  reduced  fraae  rates 
prcduce  saaller  KMSE  values  for  the  individual  fraaes  coded.  This 
indicates  that  reductions  experienced  in  f raae-to-f raae  cprrelaticns 
due  to  tenporal  subsaapling  are  completely  compensated  for  by  the 
increased  nuaber  of  bits  available  for  coding.  However,  sub jectively, 
reduced  fraae  rates  tend  to  result  in  jerky  subject  notion.  This  is 
most  apparent  for  rapidly  acving  objects  in  the  field  of  view  and  is 
of  lesser  consequence  for  slowly  changing  scenes. 

Conclusions:  Based  cn  theoretical  and  experimental  results 
obtained  to  date,  two  main  conclusions  have  been  reached.  The  first 
is  that  exploitation  of  temporal  correlations  in  addition  to  spatial 
correlations  has  been  demonstrated  to  be  a viable  technique  for  coding 
sequences  of  digital  images.  This  fact  is  demonstrated  by  a 
comparison  of  the  average  bit  rates  required  for  the  interframe 
cosine/cosine/CFCH  and  the  existing  intrafrarae  cosine/DPCM  coders  to 
achieve  the  same  level  of  NNSE  performance.  The  sixteenth  frame  of 
the  test  data  base  was  chosen  for  comparison  and  was  coded  at  an 
average  0.25  tits/pixel  by  the  interframe  cosine/cosine/DPCH  coder. 
When  using  the  intraframe  cosine/DPCM  coder,  it  was  necessary  to  code 
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this  fraae  at  a bit  rats  of  aore  than  2 bits/pixel  to  achieve  the  saae 

MUSE. 

The  second  conclusion  is  that  the  perfcraance  of  the  hybrid 
interfraae  coders  investigated  are  heavily  dependent  upon  the  type  of 
notion.  In  the  case  of  the  16  frane  head  and  shoulders  test  data 
base,  good  coding  perforaance  was  achieved  since  subject  aoveaent  was 
restricted  tc  a relatively  saall  portion  of  the  inage.  However, 
coding  perforaance  with  a different  aerial  data  base  was  degraded  froa 
the  previous  case  due  principally  to  caitera  platform  notion  which 
caused  f rane-tc-f raae  pixel  aaplitude  variations  across  the  entire 
iaage.  Since  the  perforaance  of  the  hybrid  interfrane  coders  is 
dependent  on  teaparal  correlation,  a reduced  level  of  perforaance  is 
to  be  anticipated  for  iaage  sequences  distorted  by  notion. 
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4.  Image  Restoration  and  Enhancement  Projects 

Image  restoration  and  image  enhancement  are  two  classifications 
of  image  i mj lovement  methods.  Image  restoration  techniques  seek  to 
reconstruct  or  recreate  an  image  to  the  form  it  would  have  had  if  it 
had  not  been  degraded  by  some  physical  imaging  system.  Image 
enhancement  techniques  have  two  major  purposes:  improvement  in  the 
visual  quality  of  a picture  to  a human  viewer:  and  manipulation  of  a 
picture  for  mote  efficient  processing  and  data  extraction  by  a 
machine.  Research  in  both  areas  during  the  past  si*  months  is 
described  below. 

4.1  Eigenvectors  of  Space- Variant  Point  Spread  Function  Systems 
Harry  C.  Andrews 

In  image  restoration  systems  a linear  model  results  in  an  object 
t being  mapped  into  an  imaga  ^ by  a point  spread  function  matrix  J1. 
Thus  with  noise 


Ml  + 5 L + B 


(1) 


The  simplest  linear  models  for  imaging  systems  are  given 
invariant  point  spread  functions  (SIPSF)  in  which  case 
circulant.  If  the  linear  model  is  not  space  invariant, 
represents  a space  variant  point  spread  function  (SVPSF) . In 
of  separable  systems  eg.  (1)  becomes 


by  space 
is  block 
then 
the  case 
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(2) 


G = A F B + N 

where  _A  represents  the  column  blur  and  _B  represents  row  blur  on  the 
object  .£.  In  the  SIPSF  case  ^ and  B_  are  circulants,  but  for  the  SVPSF 
case  A^  and  _E  may  have  very  general  structure.  It  is  interesting  to 
investigate  the  eigenvectors  of  such  systems  to  get  a better  feel  for 
the  underlying  eigenspace  of  the  distortions  representing  such 
systems.  In  the  case  of  SIPSF  systems,  the  eigenvectors  are  sine  and 
ccsine  waveforms  and  the  eigenspace  of  such  distortions  are  given  by 
the  Fourier  transform.  In  the  SVPSF  situation,  the  eigenvectors  often 
turn  out  to  be  variations  on  sines  and  cosines  depending  on  how 
variant  the  blur  actually  is. 

To  illustrate  thj.s  point  a separable  (SVPSF)  system  has  been 
simulated  for  two  degrees  of  blur  (moderate  and  severe) . Figure  1 
illustrates  this  situation  in  which  16  point  sources  experience 
spatially  vatiant  degradations.  The  imaging  systems  are  separable  and 
are  in  better  focus  in  the  center  and  get  more  blurred  toward  the 
edges.  Figures  2 and  3 present  selected  eigenvectors  fci  both  the 
moderate  and  severe  distortion  cases.  As  the  eigenvector  index 
increases,  the  eigenvectors  experience  an  increasing  number  of  zero 
crossings  similar  to  sine  and  cosine  functions.  Also  note  that  the 
first  eigenvector  has  no  zero  crossings  and  is  not  a constant.  These 
SVPSF  eigenvectors  appear  to  ba  FH  modulated  trigonometric  waveforms. 
It  is  interesting  to  conjecture  that  as  a function  of  the  decreasing 
variant  nature  of  the  blur  involved,  these  eigenvectors  will  converge 
to  unmodulated  trigonometric  functions.  In  examining  figures  2 and  3, 
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d)  Moderate  Distortion  Perspective 


it  is  interesting  to  nota  the  effort  each  eigenvector  goes  to  in  order 
to  resolve  finer  detail  at  certain  points  along  the  axis  compared  to 
other  positions.  Also  note  the  eigenvectors  effectively  gc  to  zetc  at 
higher  indices  in  the  center  of  the  axis  indicating  they  have  no 
effect  on  the  restoration  fcere. 


h.2  Least  Squares  Restoration  for  the  Continuous-Discrete  Model 


Steve  Hou 


For  image  restoration  purposes,  a realistic  model  is  that  given 
by  the  continuous-discrete  model  defined  by 


£= 


TI)f(e,TI)ded71 


where  a discrete  image  _g  is  obtained  from  a possibly  space  variant 

imaging  system,  described  by  hfe,^),  observing  a continuous  object 

f(e,T)).  In  digitally  restoring  such  a model  only  a finite  number  of 

* 

samples  ar«  available  for  description  of  the  estimate  f (e  ,T) ) of  the 
object.  Using  cubic  spline  interpolation 


=EECiiSi|C|Sil" 


where  Sj  (X)  is  the  i th  cubic  spline  centered  at  e. 


An  objective 


function  for  restoration  with  a smoothness  constraint  is  given  by 
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00 


W(f)  = 


IU4!i+vjJCf’'(e,Tl)]2decni 


(d) 


where 


h(e,Tl)f  (e,Tl)dedn 


(<*) 


By  lifferentiati.cn  and  subsequent  manipulation,  the  systemization 
equation  result  is 

fpTP+y  Bl  c = PT^  t5) 


Here 


P = |jH 


h (e,  "PJs  (e,  "0)  ded  Tl 


(6a) 


_sT(e,T)  = ,sT(e)®  .sV)  (6b) 

C = C.(x)  C.  (6c) 

- -iw  -j 

00 

B = j"  £n  (e,T')s',T(e,'ri)d  edn  (6d) 

-00 

Equation  t1')  is  known  as  the  normal  equation. 

The  method  of  conjugate  gradients  has  been  used  to  iteratiwely 
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search  for  the  solution  in  eq.  (5) . Because  of  computer  limitations,  a 
separable  pcint  spread  function  has  been  assumed,  for  both  space 
variant  and  invariant  systems.  For  the  separable  formulation,  the 
normal  equation  becomes 


[ A ^ A+  y B ® B 1 C=ATg 

— — — x — y — 


(7) 


where 


A =JJ  hyle, 


T)s  ( e)s  (H)  de  dn 


(«a) 


and 


/' 


^X=  Mk  <*)££  <e>dG 


(Hb) 


CO 

Vj* 

•l.  m 


V (H)  s"  (T)dli 
■i  — n 


(Be) 


f emulation,  the  generalized  extrapolated  Jacobi  iterative 
1 s-  j iwen  by 


- 1 


x ^ — y J L- 


A £-A  a+y(b  ®b  )C( 
x y 


(*) 


where  Bx  is  defined  as  Bx  except  that  no  derivatives  are  taken  of  the 
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spline  functions  in  eqs.  (8b)  and  (8c) . The  advantage  of  the 
formulation  in  eg.  (9)  is  that  no  large  matrix  inverses  need  be  taken. 

A conjugate  gradient  algorithm  has  been  icpleaented  fcr  both 
space  variant  and  invariant  cases.  The  blur  impulse  response  is  given 


H..(e,H)=hi(e)hjni) 


(10a) 


where 


h.(e)=kexp  

l o. 


( e -x  ) 


(10b) 


and  a.  =|kx.|  such  that  k governs  the  amount  of  blur  or  spread  as  a 

function  of  position  (x^)  in  ^e  imaging  plane.  A similar  equation 

results  for  h.  (T>) . For  the  space  invariant  case  a.  was  set  equal  to  k 

without  x.  contributing  to  the  spread  of  a.  . 

x 1 

The  simulated  results  by  using  the  conjugate  gradient  algorithm 
are  shown  in  figures  1 through  6.  For  both  restorations  from  moderate 
sipsf  blur  (figure  1)  and  from  moderate  SVPSF  blur  (figure  2),  the 
results  are  strikingly  good  for  Y =0.  The  justification  foe  such 
results  is  that  the  PSF  is  fairly  localized  (i.e.  narrow) , and  thus, 
the  matrix  fli  is  well  conditioned.  In  other  words,  the  eigenvalues  of 
A.  are  clustered  together  sc  that  A.  is  far  from  singularity. 


On  the  other  hand,  as  the  PSF  spreads  out  and  the  image  becomes 
more  blurred,  the  restored  objects  for  both  SIPSF  and  SVPSF  are  far 
from  perfect.  For  x=0  ringing  in  separable  form  shows  up  in  the  SIPSF 
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(a)  Original 


Figure  4.  2-1.  Restoration  from  moderate  SIPSF  blur  (k  = 1) 


rdr*- 


(b)  Restored  C for  5 = 0 


Figure  4.  2 


(c)  Restored  F for  6=0 


6 . Restoration  from  very  severe  SVPSF  blur  (k  = 1) 


case;  and  the  norm  of  the  error  matrix  in  the  gradient  algoritha 
oscillates.  This  is  because  the  matrix  Jk  now  is  badly  conditioned  and 
approaches  singularity.  Theoretically*  as  Y=0,  the  conjugate  gradient 


P 

* 

H 

:'S 

f 

% 

f 


algoritha  is  the  saae  as  the  pseudoinverse  of  _A.  Under  this 
condition*  the  ellipsoidal  contour  surface  in  the  direction 
corresponding  to  zero  eigenvalues  shrinks,  thus  residual  errors  oan  no 
longer  maintain  orthogonality,  and  the  computing  tiae  to  convergence 
gccws  enormously. 

As  shown  in  figures  1 and  2,  the  tradeoff  between  the  picture 

smoothness  and  sharpness  which  may  be  accompanied  by  oscillations 

-4 

becomes  evident  frca  the  results  for  Y =10  and  Y =0  in  both  SIESF  and 

SVPSF  cases.  The  price  paid  for  sharp  pictures  is  a long  iteration 

-8 

time.  Notice  that  in  the  SIPSF  case,  the  restored  object  for Y =10  is 

almost  identical  with  that  for  Y=0.  Hence,  it  is  suspected  that  in 

-6 

the  SVPSF  case,  the  oscillation  could  be  supressed  by  using  Y = 10  or 
-7 

10  without  much  impairaent  of  the  picture  sharpness,  but  with  the 
additional  advantage  of  faster  convergence. 

A 

Tha  white  spcts  appearing  in  all  the  C pictures  are  the  negative 

A 

coefficients  in  the  C matrix.  Because  of  the  positive  nature  of  the 
spline  basis  function  the  coefficients  must  have  negative  values  in 

A 

order  to  reconstruct  f(e,Tl)  properly.  As  expected,  the  white  spots 
appear  at  the  high  contrast  areas  of  the  GIRL  picture,  such  as  along 
the  edges  of  her  scarf,  or  the  flowers  and  in  her  eyes.  As  Y 
decreases,  the  number  of  white  spots  increases  because  the  restored 
picture  becomes  sharper.  Foe  severely  blurted  images,  the  white  spots 


are  scarce  and  hence,  the  object  is  no  logger  sharply  reconstructed. 

4.3  A General  Image  Estimation  Algorithm  Applicable  to  Multiplicative 
and  Non-Gaussian  Noise 

Nasser  E.  Nahi  and  Kohanned  Naraghi 

In  statistical  iaage  enhancement,  an  image  is  described  by  a 
two-dimensional  random  process  (field).  These  processes  are  often 
characterized  by  their  mean  and  autocorrelation  [3],  Denoting  the 
image  brightess  function  by  b{i,J),  with  i and  j as  the  horizontal  and 
vertical  variables,  the  twc  moments  are  defined  as 

M(i,j)  = E{b(i,j)}  (1, 

R(i>  j.  k.i ) = E{  fb(i,  j))-M(i,  j)][b(k, /)-M(k,  t )]1  (2) 

where  E is  the  mathematical  expectation  operator.  The  degraded  image 
(ccmmonly  referred  to  as  the  observation)  is  denoted  by  y(i,j)  and 
specifies  the  functional  relationship  between  signal,  h(i,j),  and 
noise  (i,j)  given  fcy 

y(i,  j)=  f [b(i,j),Y(i,  j)}  (3) 

where  f may  be  ncnlinear  andY  (i,j)  may  be  vector  valued. 

Optimum  filtering  of  images  under  the  general  condition  of  eg.  (3) 
has  received  little  attention.  However,  a variety  of  procedures  have 
been  developed  for  the  special  linear  case,  where 


y(i»  j)  = b(i,  j)+  Y(i*  j) 


<«*) 


* 


where  Y (i,  j)  is  white  and  Gaussian  [11  to  17],  Although,  eg.  (4) 
describes  many  natural  forms  of  degradations  [12  to  16],  there  are 
conceivably  as  nany  situations  where  this  model  does  not  apply. 
Examples  are  iiages  with  film  grain  noise  and  pictures  cfcserved 
through  non-hcmogenepus  clcud  layers,  where  the  noise  is  a random 
attenuation  factor.  In  these  examples,  the  observations  take  the  form 

y(i»  j)  = Y (i.  j)b  (i,j)  (5) 

The  majority  of  the  existing  linear  estimation  procedures  require 
the  correlation  function  3(i,j,k,l)  to  be  specified  as  an  analytic 
function  of  a particular  form  [12  to  17].  This  limits  the  generality 
of  these  methods  since  they  cannot  be  applied  to  practical  cases 
where,  the  function  R{i,j,k,i)  is  often  specified  numerically  at  only 
a small  number  of  argument  indices. 

The  purpose  of  this  week,  is  to  develop  a general  estimation 
method  which  requires  numerical  values  of  the  autocorrelation  function 
only#  and  is  applicable  to  nonlinear  (as  well  as  linear) 
observation  systems.  Furthermore,  the  estimation  technique  will  be  of 
recursive  nature,  and  hence,  computationally  efficient. 

Notation;  An  image  is  viewed  as  an  n x n matrix  with  elements 
b(i,j),  where  b(i,j)  is  the  intensity  cf  the  image  at  pixel  (i,j).  To 
reduce  notaticnal  complexity  the  pixels  are  indexed  by  1,2,...,n 
consecutively  ficm  left  to  right  and  top  to  bottom.  This  convention 
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enables  reference  to  the  doubly  indexed  image  b(i,j)  as  b(k), 
symbolically.  Hence  egs.  (1)  to  (1)  can  be  written  as 


M(k)  = E[b(k>}  (6) 

R<M)=  E{fb(k)-  M(k)]fb(/)-  MW]]  (7) 

y(k)  = f{b(k),  y (k)]  (3) 


Let  the  process  *(k)  be  defined  as 

x(k)=  b(k)-M(k)  (9) 


for  k-1,2,...,n  . Thus,  the  problem  of  estimating  b(k)  reduces  to 
estimating  x (k) . 

Estimation  Method:  The  minimum  mean  square  (MHS)  estimate  *a(k), 
of  a process  x(i)  at  time  (pixel)  k and  for  a given  set  of  observation 
yfl)  (k)  is  given  by  £233 

xG(k)=  E [x(k)  | y(l y(k)}  (10) 


Letting 


Y(k)=[y(l),...,y(k-1)}  (11) 

2 

then  it  can  be  shown  [21,25]  that  x (k)  and  its  error  variance  o°  (k) 
are  functionally  related  to  7(k)  and  y (k)  by 
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xa(k)  = 


jx(k)p(y(k)  | x(k))p(x(k)  | Y (k))dx(k) 
Jp(y(k)  |x(k))p(x(k)  I Y(k))dx(k) 


(12) 


.2  [[x(k)-xa(k)l  p(y(k)  |x(k))p(x(k)  | Y(k))dx(k) 

- -s/L J 


0°  (k) 


j”p(y(k)  |x(k))p(x(k)  | Y(k))dx(k) 


(13) 


where  p designates  appropriate  density  functions.  Equations  (12)  and 
(13),  in  turn  suggest  that  the  optimal  estimation  at  k is  achieved  by 

first  finding  p(x(k)t?(k))  and  then  using  it  along  with  y { k y to  arrive 

2 

a a , 

at  x (k)  and  a (k)  . The  mean  of  p(x(k)|Y(k))  is  the  MMS  one  step 

prediction  of  the  random  variable  x(k)  and  its  ’'ariance  is  the  error 

variance  of  the  predicted  value.  Thus,  the  optimal  estimation  at  time 

k can  be  thought  of  as  a two  step  procedure  depicted  in  figure  la, 

where  blocks  P and  F may  be  identified  as  the  prediction  and  filtering 

steps,  respectively.  In  this  system  structure,  y (k)  is  isolated  from 

other  random  variables  and,  assuming  p(x(k)|t(k))  is  Itncwn, 

conceptually  one  can  deal  with  its  ncnlinearities  in  block  F,  i.e.  if 

2 

p (x  (k)  I Y (k) ) is  given,  then  derivation  of  x°(k)  and  oa  (k)  is 
accomplished  by  carrying  out  the  integrations  in  egs.  (12)  and  (13). 
However,  for  the  general  observation  of  eg.  (3),  derivation  of  this 
probability  density  does  not  lend  itself  to  analytic  methods  and 
available  numerical  approaches  are  computationally  unfeasible  [23, 
Chapter  7 ]. 

In  this  report  an  alternate  procedure  is  considered,  whereby  an 
approximation  to  the  protatility  density  p(x(k)|Y(k))  is  derived*  The 
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method  is  compatible  with  the  logic  of  the  estimator  in  figure  la. 
This  logic  consists  cf  representing  past  information  (i.e. 
information  due  to  a priori  statistic  and  observations  y (1 ) , . . . , f (K- 1) 
in  the  form  of  a probability  density  to  be  combined  with  y ( k)  in  block 
F.  Based  on  this  premise  and  the  goal  of  algorithmic 
implementability , the  estimator  is  constructed  according  to  the 
following  restrictions. 

a) .  Only  the  first  two  moments  of  any  random  variable  are 
computed. 

b)  . The  prediction  process  is  chosen  to  be  linear. 

c)  . The  prediction  is  to  be  based  on  a selected  small  number  of 
past  estimates.  This  will  impose  a desired  limited  memory 
requirement  for  the  estimator. 

A A 2 

Letting  x (i)  and o (i)  represent  the  estimate  and  its  error 
variance,  respectively,  at  time  i,  then  the  block  diagram  in  figmre  1b 
represents  the  structure  of  the  proposed  estimator.  In  this  figure 
blocks  LP,  F and  C signify  linear  prediction,  filtering  and  cne  unit 
time  delay,  respectively.  The  subscript  R is  an  indication  of  the 
size  of  memory  and  x':{h)  and  a (k)  are  the  one-step  predicted  value 
and  its  error  variance.  The  set  { k-I j,. . . ,k-lM]  is  a set  of  two 
dimensional  indices  each  distinct  and  prior  to  k. 

Modeling  Procedtre;  To  derive  the  linear  predictor  (block  IP  of 
figure  1b),  the  a priori  correlation  information  is  first  incorporated 
into  a linear  finite  order  model  of  the  process  x(k)  in  the  form  of 
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(14) 


x(k)  =^Y3ix(k“1i 


) + B u(k) 


where  0,  constants  and  {u  (k)  ,u  (k—  1 ) ,.  ..  1 is  a set  of 

independent  identically  distributed  random  variates  with 


E{u(k)]  = 0 
E[u(m)u(n)}  = 


0 if  m i n 

1 if  m = n 


(15) 


Consequently,  e ij . . { 1 4 ) is  an  autoregressive  model  [12],  [18  to  20]. 

The  problem  of  modeling  consists  of  determining  the  order  PI*  the 


coefficients  Pj  ,...,£ 


M' 


the 


set 


of  two  dimensional  indices 


k-Ij , . . . , k-IM  and  the  variance  of  the  white  noise  term  B u (k)  in  eg. 
(14).  In  this  work,  first  a procedure  is  developed  to  derive  an 
autoregressive  model  for  a given  M followed  by  a discussion  cn  the 
best  choice  of  P.  The  modeling  criterion  is  chosen  to  be  minimization 
of  e [ B u(k)}.  The  procedure  uses  the  numerical  values  of  the 
correlation  function  and  does  not  require  analytic  representation  of 
R(m,n).  The  results  are  illustrated  by  the  following  example. 

Consider  the  stationary  two-dimensional  correlation  function 


R(i,  j,k,2)  = R(|i-k|,  |j-2  |)=  E{x(i,  j)  x(k,2)}  = exp  £ - J (i-k)2  + (j-2)2  J 


Application  cf  above  procedure  provides  the  following: 

a)  . Best  2nd  order  modal  is 

x(i,  j)  = 0.  3 x(i,  j-1)  + 0.  3 x(i-l,  j)  + 0.  883  u(I,  j) 

b)  . Best  3rd  order  mcdel  is 

x(i , j ) = 0.29  x(i,  j-1)  + 0.25  x(i- 1 , j)  + 0.12  x(i- 1,  j+1)  +0.87?  u(i,  j) 
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c) . Best  4th  order  model  is 


x(i,  j)  = 0.  28  x (i,  j-1)  + 0.  24  x(i-l,  j)  +0.  12  x(i-l,  j+1) 

+ 0.03  x(i-l,  j-1)  + 0.8769  u(i,j) 

d)  . Best  *th  order  model  is 

x(i,j)  = 0.28  x(i,  j-1)  + 0.24  x(i-l,  j)  +0.11  x(i-l,  j+1) 

+ 0.03  x(i-l,j-l)  + 0.02  x{i=l,j+2)  + 0.8768  u(i,  j) 

Hence,  for  example,  to  a third  decimal  place  accuracy,  the  3rd  order 
model  is  a sufficient  approximation.  Note  that,  for  example,  the 

derivation  of  the  3rd  order  model  requires  the  numerical  values  of 
R(C,C),  R (0 , 1 ) , P ( 1 , 0)  and  R(1,1). 

Linear  Prediction:  Let  the  model  of  the  random  process  x (It) 
(obtained  in  the  previous  section)  be 

M 

x(k)^8.  *(k-I.)  + Bu(k)  (16) 

i=l 

A 

Given  the  estimate  x (i) , i = 1 , 2, . . . , k- 1 the  linear  prediction  x (k) , in 
general,  is  given  ty 

k-1 

x' (k)  aj  x(k-j)  (17) 

j = l 

vhereCLj  ,...,a^  ^are  to  be  chosen  such  that 

E[ x(k)  - x::  (k)l2  (18) 
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is  minimized.  This  minimization  is  to  be  carried  out  subject  to  the 
system  structure  of  figure  1b  and  is  based  on  available  infpraation  to 
the  predictor.  This  inforaation  consists  of  the  values  of  x (i)  and 

A A 2 

(i)  , i<k-1.  Since  each  x ( 1 ) anda  (i)  is  the  *ean  and  variance, 
respectively,  of  a posterior  density  on  x (i)  at  tiae  i (having  used 
observations  through  y (i)) , then  the  expectation  in  eg.  (13)  is  sell 
defined  and  operates  on  each  randoa  variable  x(i)  such  that 


E{x(i)}  = x(i) 


E{[x(i)  - x(i)]2]  = o2(i) 


(19) 


Theorem  1:  Shen  the  random  process  x(k)  satisfies  eg.  (16) , then  the 
(optimal)  choice  of  ctj  , ,. . . ,a^ jin  eg.  (17)  which  minimizes  eg.  (19) 
is  given  by 


a. 


^ if  k-j  = k-L 
0 otherwise 


The  proof  is  given  in  [26]. 

This  theorem  states  that  the  best  linear  predictor  is  given  as 

M 

x"(k)  =y^B.x(k-l.)  (20) 

fel 
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Thi 


The  implexe ntation  of  eq.  (20)  is  very  simple.  This  simplicity, 
along  with  the  effectiveness  of  the  result  as  illustrated  in  the  next 
sections,  ace  the  justification  behind  the  necessary  approximations. 

Filtering  Step:  Referring  to  figure  1b,  the  computat icnal  logic 

of  block  F is  now  developed.  The  predicted  value  x (k)  and  its 
2 

variance  0^  (k)  , obtained  from  the  linear  predictor,  represent  the  mean 
and  variance  of  the  a posteriori  density  on  x(k).  This  density 
represents  the  available  knowledge  on  the  random  variable  x(k)  prior 
to  reception  of  y(k).  Since,  for  a given  mean  and  variance  the  normal 
distribution  represents  the  maximum  uncertainty  (entropy)  [24.,  p. 

132],  this  density  function  is  assumed  to  be  normal.  Further 

*2 

uncertainty  is  associated  kith  x (k)  if  o (k)  is  used  in  place  of 

2 

0^  (k) . Consequently,  an  approximate  and  a rather  conservative  choice 
of  the  probability  density  for  x(k)  is 


Pr*(k)l  = fa'(k)  J 2 tt]"  exp|- 


rxoo-xV)]2 


2c*2(k)  ) 


(21) 


Observation  y(k)  and  p(x(k))  in  eg.  (21)  are  combined  to  derive 

A 

the  Bayes  estimate,  x(k) 


x(k)  = E{x(k)  jy(k)}=  J x(k)  p(x(k)(y(k))  dx(k) 
= p(y(k^  |x(k)  p(y(k)  | x(k))p(x(k))  dx(k) 


(22) 


I 

i 


■aagma; 


But 


f 


t 


if 

* 

% 

% ■ 

; T' 

*• 

4; 

% 

.♦ 


p(y(k))  = j P(x(k),  y(k))  dx(k)  = J" p(y(k)  | x(k))  p(x(k))  dx(k) 


Hence,  [ 25  ] 


x(k)  = 


J 


x(k)  p(y(k)  \ x(k))  p(x(k))  dx(k) 


f 


p(y(k)  | x(k))  p(x(k))  dx(k) 


(23) 


Similarly, 


S2(k)  = E|fx(k)  -x(k)]2  |y(k)|  = 


j(x(k)-x(k))2p(y(k)  |x(k))p(x(k))dx(k) 


jp(y(k)l 


x(k))p(x(k))  dx(k) 


(24) 


where  p(x(k))  in  eqe.  (23)  and  (24)  are  given  by  eg.  (21)  and 
p (y  (k) | x (k) ) is  obtained  from  the  observation  system  structure. 


a 2 

In  general,  evaluation  of  x (.)  and  o (.)  in  eqs.  (23)  and  (24) 
will  te  performed  numerically.  This  in  turn,  allows  the  procedure  to 
be  applicable  to  a broad  class  of  observation  systems  including 
nonlinear  forms  cf  the  observation  y ( k) . The  feasibility  of  this 
estimator  is  due  to  the  structure  of  figure  1b  which  leads  to  egs. 
(23)  and  (24)  . 


Multiplicative  Noise  Term  in  Observation;  Consider  observations 
containing  uniform  multiplicative  noise.  In  this  case  the  observation 
is  given  by 

y(k)  = Y(k)  fx(k)  •(-  M(k)]  (25) 


1 


i 


i 


with 
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if  0 < Yj(k)  * Y(k)  < Y2( 


P(Y(k))  = 


otherwise 


With  x (k)  ♦ n (k)  as  the  image  intensity  at  pixel  k,  eqs.  (23)  and 
(24)  become  [ 25  ] 


v 1 C x(k) 
x(k)  = G ) x(k)+M(k)  eXP 

a 


(x(k)  - x"'(k))‘ 

2c*2(k) 


dx(k) 


a 2 j f fx(k)-x(k)]  (x(k)-x' 

0 (k)  =G  3 x(W  * M(k)  e!£P  - ~Za'nm 


(x(k)  - X (k))“1 

exp r? |dx(k)  (28) 


where 


1 x(k)  -x  (k)) 

exp  - rj dx(k) 

x(k)  + M(k)  2o"  (k) 


= -m.  _i 

Y2(k) 

= -likI  . 

Yj(k) 


Since  egs.  (27) tc  (29)  are  definite  integrals,  they  can  be  evaluated 
numerically.  All  noisy  images  contain  uniform  multiplicative  noise 
with  noise  bounds  as  indicated  in  these  figures.  The  estimated  images 
of  figure  2 to  4 provide  5.48,  7.50  and  7.7  db.  improvement. 


respectively.  Aside  fecit  this  guantitative  improvement,  the 


Figure  4.3-3.  Uniform  multiplicative  noise 


(b)  Noisy,  noise=0.7-l  (c)  Estimate 


Figure  4.  3-4.  Uniform  multiplicative  noise 

r 

* 

* 

A * 
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preservation  of  edgas  in  the  estimated  images  should  be  noted.  The 
responsiveness  of  the  estimator  to  abrupt  pixel  tc  pixel  intensity 
changes  is  due  to  the  estimator  structure  of  figure  1b. 

The  estimation  procedure  can  also  be  applied  to  more  general 
observation  systems.  As  an  example  consider  the  case  where 

y(k)  = y(k)  [x(k)+  M(k)]  + v(k)  (31) 

where  Y (k)  and  v (k)  are  both  uniform.  Letting  the  density  of  Y(k)  be 
given  by  eg.  (25)  and  that  cf  v(k)  be 


p(v(k))  = 


1 

v2(k)-v1(k) 

0 


if  Vl(k)  £ v(k)£  v2(k) 
otherwise 


(32) 


then  p(y(k)lx(k))  can  be  obtained  in  terms  of  the  convplutipn  of 
p(Y(k))  and  p(v(k))  [22].  This  density,  then,  can  be  substituted  in 
egs.  (23)  and  (24)  to  obtain  pertinent  filtering  eguations  [25]. 
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4.4  Image  Restoration  by  Smoothing  Spline  Functions 
Mohammad  J.  Peyrovian  and  Alexander  A.  Sawchuk 

In  a linear  space-invariant  imaging  systea  with  ppint-tpread 
function  h (x) , the  iaage  g (x)  is  given  by 

g(x)  = Jh(x-u)  f{u)du+n(x)  (1) 

where  n (x)  represents  aeasureaent  Noise.  In  order  to  estimate  the 
object  function  f (u)  froa  image  g(x)  by  a digital  computer,  the  above 
continuous  model  must  be  discretized.  A common  method  is  to  sample 
the  functions  h and  g at  a finite  number  of  points.  Spline  functions, 
because  of  their  highly  desirably  interpolating  and  approximating 
characteristics,  are  an  interesting  alternative  to  the  above  method. 
For  uniformly  spaced  knots,  a class  of  spline  functions,  called 
B-splines,  has  the  following  properties 

(i)  shift  invariance 

(ii)  strictly  positive 

(iii)  convolutional  property 

(iv)  local  basis  property 

Using  B-splines  for  interpolation  or  approximation,  the  functions 
f and  h can  be  represented  by  B-splines  of  degrees  m and  n, 
respect ively 
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12) 


f(x)  =2^  fi  Bm  {x-xi> 
i=-06 


h(x>  = Ehj  B„(x-xj)  (3> 

j=.oo 

Substituting  egs.  (2)  and  (3)  in  the  convolution  integral  of  eg.  (1) 
gives 


CO  00 


hi  B (x-x  )*B  (x-x  ) + n(x) 
3 m in  j 


From  the  convc  luticnal  property  of  B-splines 


(4) 


B (x-x.)*  B (x-x.)=  B , (x-x.-x.)  (5) 

m'  r n'  y m+n'  i y 

and  representing  g (x)  by  B-spline,s  of  degree  tn  ♦ n and  assuaing 
A x=xi+1-Xj  gives 


2gkBm+n(x-kix)=2  Bm+n(x-|i+j),'x)+n<x) 

k=  -00  i=-00  -CO 

Eguations  (4),  (5)  and  (6)  show  that  the  B-spline,  which  is 

interpolating  the  deternin istic  part  of  the  degraded  inage,  nust  be  of 
higher  degree  than  the  B-splines  interpolating  object  and  point-spread 
function.  In  ether  words,  since  the  blurred  iaage  is  always  saoothec 
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than  the  object,  a higher  degree  spline  can  follow  the  image  function 
better  than  the  one  approximating  the  object  function.  This  oan  be 
explained  in  the  Fourier  domain  by  observing  that  the  Fourier 
transform  of  an  i-th  degree  B-spline  is  a Sine  function  to  the  power 
a.  As  m increases  the  aaplitude  of  highar  frequencies  decreases. 
Since  a blurred  iwage  has  less  higher  frequency  content  than  the 
object,  a higher  B-spline  can  represent  the  inage  better  than  the  one 
representing  the  object. 

In  a noiseless  imaging  system,  eg. (6)  may  be  written  in  the 
matrix  form 


g = H f 


(7) 


If  the  point  spread  function  is  of  finite  width,  the  matrix  11  is 
banded.  Figure  la  is  a rectangular  object  which  is  blurred 
analytically  by  a 4th  order  polynomial 


h(x) 


-3.  5 Sx^3.5 


(8) 


, elsewhere 


The  object  is  a stop  function,  therefore  it  is  interpolated  by  a zero 
order  B-spline.  The  second  derivative  of  h at  points  x=-3.5  and  x=3.5 
is  a step  function  and  it  is  interpolated  by  a second  order  B-spline. 
Since  the  convclution  of  a zero  and  second  order  B-spline  is  a cubic 
B-spline,  the  image  is  interpolated  by  a cubic  B-spline.  Figure  1b, 
the  restored  image  with  and  without  splines,  shows  that  the  spline 


mu 
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Spline  restoration. 


restores  the  edges  such  sharper  than  the  coason  pulse  approximation 


■ethod.  Figure  2 is  another  exauple  of  spline  restoration  applied  to 
a two  diuensional  blur  with  point  spread  function 

H(x,y)  = h(x)h(y)  (9) 

where  h is  defined  in  eg. (9). 

For  a noisy  iuage,  the  iaage  data  is  first  saoothed  by  ainiaizing 


i 


1 


[g"(x) ]2 


d x 


aaoiig  all  functions  gee  such  that 


(10) 


Here  y-  is  the  noisy  iaage  neasured  at  point  Xj...s>0  and  o.  >o  are 
given  nuabers.  Setting  S-0  leads  to  an  interpolation  problem.  The 
factor  o.  control  the  saoothing  window  at  point  x.  and  S controls  the 
extent  of  saoothing.  If  the  standard  deviation  of  y.  is  available,  it 
aay  be  used  as  a . In  this  case,  natural  values  of  S lie  vithin  the 
confidence  interval  of  the  left  hand  side  of  eg.  (10)  as  given  by 


i l 

N - (2N) 2 SSSN  + (2N)2 


where  N is  the  nuaber  of  data  points.  Reinsch  [3]  has  shown  that  the 
solution  to  egs.  (9)  and  (1C)  is  a cubic  spline,  and  more  generally,  is 
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(a)  Original  image 


■ U5E 


(b)  Blurred  image  (c)  Restored  image  using 

spline  functions 


Figure  4.4-2.  Examples  of  spline  restoration 
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a spline  function  of  degree  2K- 1 for  least  sguare  minimization  of  the 


K-th  derivative  instead  of  the  second  derivative.  In  smoothing  (S>0) , 
the  shape  of  the  function  is  such,  sore  influenced  by  the  minimum 
principle  of  eq.  (9)  than  in  interpolation  (S=0) . 

The  above  saoothing  criterion  will  be  subject  of  further  research 
on  noisy  blurred  inages,  particularly  the  case  K=2  because  it  leads  to 
cubic  splines  which  read  sinpler  algorithms  and  less  coaputation. 
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4.5  Detection  and  Estimation  of  Iaage  Degradedby  Film-Grain  Noise 


Firouz  Naderi  and  Alexander  A.  Sawchuk 


The  goal  cf  this  research  has  been  to  analyze  the  profalea  of 
fila-grain  noise  in  the  context  of  detection  and  estiaation  theory. 
The  first  step  is  the  development  of  a mathematical  acdel  that 
reflects  some  of  the  complexities  of  image  formation  process,  and  yet 
is  tractable  in  the  subsequent  restoration  of  the  image. 
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Denoting  by  y(i,j)  the  observed  optical  density  of  photographic 


fila  as  aeasured  by  a aictc  densitcaeter  let 


y(i.  j)  = s(i,  j)+  n(i,  j) 


(1) 


where  S(i,j)  denotes  the  density  that  Mould  have  been  registered  in 
the  absence  of  grain  noise  and  n(i,j)  is  the  noise.  Experiments  by 
researchers  in  the  field  of  Photographic  Science  have  indicated  that 
n(i,j)  is  approxiaately  Gaussian  distributed  with  zero  aean  and  a 
variance  that  is  dependent  on  the  type  of  the  filas  used,  the  siae  of 
the  scanner  apertare  and  the  value  of  S(i,j).  Clearly  the  pbservation 
aodel  described  in  eg.  (1)  is  additive  with  signal-dependent  noise. 
Equivalently,  the  additivity  of  this  aodell  aay  be  sacrificed  to 
obtain  a signal-independent  noise  aodel.  The  result  of  doing  so  is 
the  nonlinear  observation  aodel 


y(i,j)  = S(i,j)  + gfs(i,  j)]  n(i,  j) 


(2) 


where  the  noise  n.(i,j)  is  zero  aean  and  unit  variance  Gaussian.  The 
fora  of  the  function  g(.)  has  been  subject  cf  soae  discussion.  The 
experiaental  fora 


g[S(i,j)]=  k[s(i,j)]b 


(3) 


Jj 


\ 

r 
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has  been  found  to  be  in  agreement  vith  aany  different  theoretical  and 
experiaental  results.  Simplified  photographic  emulsion  models  such  as 
Nottings,  result  in  a value  of  1/2  for  the  exponent  b in  the  above 
equation.  Data  taken  by  Higgens  and  Stultz  [ 1 ] suggest  values  of  b in 
the  range  0.3  to  C.4  if  the  scanning  aperture  is  allowed  to  vary 
within  a reasonable  range. 

With  this  model  the  restoration  problem  is  considered  in  two 
different  contexts:  detection  and  estimation.  In  aany  image 

processing  problems,  it  is  necessary  to  use  a high  magnification  to 
extract  image  information  out  of  a photographic  recording.  A digital 
image  of  size  256  x 256  can  be  obtained  by  scanning  a square  region  of 
side  approximately  1.25  mm  using  a 5 micron  aperture.  Measuring 
optical  density  in  such  a small  region  of  a photographic  film  results 
in  such  a high  level  of  grain  noise  that  distinguishing  between 
adjacent  areas  of  small  contrast  with  the  naked  eye  becomes 
impossible.  Recently  Zueng  and  Barrett  considered  image  detection  by 
a method  called  the  "Noise  cheating  algorithm."  References  [2,3]  show 
that  this  equivalent  to  method  is  sub  optimal  maximum  liklflhood  j 

detection. 

To  set  up  tbe  problem  in  the  framework  of  detection  theory 
suppose  that  the  portion  of  the  photographic  film  which  is  to  be 
scanned  can  be  segmented  into  n spatially  uniform  or  near  uniform 
density  regions  R 1 , . . . , R^p  Let  a square  aperture  of  size  a x a be 
used  to  measure  the  optical  density  of  the  fils.  It  is  then  possible 
to  formulate  an  M ♦ 1 hypothesis  problem.  The  first  K hypotheses,  Hj 
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are  the  hypotheses  that  a given  densitometer  reading  was  obtained  «hen 


E 

fi 

;i  the  aperture  was  entirely  in  one  of  the  fi  regions  R.  . The  last 

!, 

' hypothesis  HM+1  corresponds  to  a reading  taken  when  the  aperture 

overlapped  on  two  or  aore  regions  siaultaneously  as  shown  in  figure  1. 
Conventional  aaxiaui  likelihood  or  Bayesian  detectors  can  now  be 
utilized  for  optiaal  detection  of  the  M ♦ 1 hypotheses. 

A siaple  suboptiaal  aethod  to  accoaplish  this  procedure  is  to 
perfora  the  M ♦ 1 hypothesis  detection  in  two  different  steps.  In 

I step  one  the  hypotheses  ^M+1is  ignored  and  the  other  H hypothesis  are 

t 

I optimally  detected.  Therefore,  in  the  first  step  the  possibility  that 

some  readings  aight  have  been  taken  when  the  aperture  overlapped  aore 
than  one  region  is  ignored.  In  the  second  step , in  regions  when 

hypothesis  HM+1  appears  to  be  highly  probable  (i.e.  the  edges),  the 
iaage  is  re-examined  with  a finer  aperture  to  recover  details. 
Figures  2c  to  2e  contain  simulation  results  of  this  restoration 
procedure  for  the  three  detection  strategies  described  below. 

ttaximua  lihlihcod  detection  Ecr  signal-independent  noise: 

Referring  to  figure  1 assume  that  the  aean  density  in  region  R , 

called  the  background,  is  and  the  variance  of  the  readings  taken 

2 

with  an  aperture  of  size  a x a in  this  region  is  G^.  The  scanned 
image  is  of  size  256  * 256.  A two  by  two  spatial  averaging  is  first 
performed  on  the  scanned  image  (Note  that  in  effect  the  averaged  image 
is  what  we  would  have  obtained  had  we  scanned  the  film  with  a 2a  x 2a 
apreture  to  begin  with.)  In  the  averaged  image,  pixels  in  the  region 
I Rj  will  now  have  mean  m^  and  variance  G^/4.  Each  pixel  in  the 

I 
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(b)  Image  with  film- grain 
noise  added 


(d)  Maximum  likelihood  detected 
image  assuming  signal 
dependent  noise 


(c)  Maximum  likelihood  detected  image 
assuming  signal  independent  noise 


(e)  Bayesian  detected  image 


Figure  4.  5-2  , Image  detection  in  the  presence  of  film-grain  noise 


i 


\ 


averaged  image  is  new  quactized  to  one  of  'I  levels.  These  levels  are 
chosen  such  that  one  of  them  will  coincide  with  the  mean  of  the 
background,  m^#  and  the  others  will  be  4 apart  from  each  pther. 
Since  the  distribution  of  the  noise  is  Gaussian,  if  th*  decision 
levels  for  the  quantization  are  set  exactly  at  the  aid-point  between 
each  quantization  level  then  it  is  easy  to  deaonstrate  that  the 
quantization  is  in  fact  maximum  liklihood  detection. 

Since  the  levels  are  taken  to  be  four  standard  deviations  apart, 
all  the  image  regions  which  happen  to  have  a mean  density  equal  to  one 
of  the  quantization  levels  will  almost  always  be  restored  to  their 
correct  mean  density  following  the  quantizaton.  Regions  having  mean 
densities  that  fall  between  two  quantization  levels  will  te  "coded" 
into  a percentage  of  these  two  levels. 

The  second  step  in  the  maximum  liklihood  detection  process  is  to 
rework  the  edges  in  the  quantized  image  by  comparing  the  quantized 
image  with  the  original  scanned  image  which  was  scanned  with  the  finer 
a x a aperture.  Figure  2c  is  the  detected  image  using  this  procedure. 

Kaximun:  liklihood  detection  for  signal-dependent  noise:  The 
performance  of  the  previous  detector  is  dependent  upon  the  distance 
between  the  quantization  levels.  If  the  levels  are  four  standard 
deviations  apart,  it  is  certain  that  regions  whose  mean  densities 
coincide  with  cne  of  the  quantization  level  will  be  clear  of  the 
noise.  As  seen  in  eg.  3,  the  standard  deviation  of  the  noise  is  a 
function  of  the  signal.  Therefore  for  an  image  with  high  dynamic 
range  it  is  necessary  to  increase  the  distance  between  the  higher 
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quantization  level  so  as  to  keep  the  listance  always  four  G . 
Furthermore,  since  the  standard  deviation  varies,  the  decision  level 
of  the  quantization  which  corresponds  to  the  maximum  likelihood 
detection,  will  no  lpnqer  te  at  the  aid-point  between  the  qmantization 
levels.  Figure  2 shows  the  iuproveaent  over  the  previous  detector 
when  the  signal  dependence  of  the  noise  is  taken  into  account  with  the 
proper  quantization. 

Bayesian  Cetection:  As  previously  mentioned,  quantization  is,  in 
effect,  maximum  likelihood  detection.  To  take  advantage  of  any  a 
priori  knowledge  that  might  be  available  about  the  image,  it  is 
advantageous  to  perform  Bayesian  detection.  Corresponding  to  the  N 
hypothesis  detection  in  the  first  step  of  the  above  two  detectors,  the 
mean  densities  of  the  fl  region  may  assume  a distribution  over  a small 
range.  Using  the  distribution  as  apriori  statistics,  the  result  of 
fcayesian  detection  is  shown  in  figure  2e. 

Summary:  Estimation  algorithms  are  presently  being  applied  to 
film-grain  noise.  Both  Wiener  filter  and  a nonlinear  filtering 
reported  in  BSC  image  processing  institute  report  580  [4],  will  be 
applied,  and  their  performances  will  be  compared  and  reported. 
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4.6  Vignetting  and  Density  Correction  for  CRT  Film  Recording 


Herner  Frei 


The  acquisition  of  digitized  image  data  and  the  restitution  of 
processed  pictures  are  generally  costly,  time-consuming,  and  yet 
essential  steps  of  digital  image  processing.  Errors  and 
non-linearities  introduced  by  the  scanning  and  display  equipment  or 
the  photographic  process  can  add  a surprising  amount  of  unwanted  and 
uncontrolled  "image  processing."  These  parasitic  effects  are  by  no 
means  always  readily  visible  in  the  finished  product,  but  they  may 
well  invalidate  the  results  of  computer  image  manipulations.  a 
careful  control  of  the  electro-optical  machinery,  the  photographic 
process,  as  well  as  an  understanding  of  human  visual  factors  is 
therefore  essential  to  inscre  the  success  and  credibility  of  digital 
image  processing. 
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Visual  Factors;  Optimum  reflection  prints,  transparencies  and 


television  images  practically  never  replicate  the  brightness 
distribution  at  original  scenes,  in  the  sense  that  color  images  do  not 
reproduce  the  spectral  energy  distribution  of  colored  lights. 
Although  comprehensive  fidelity  criteria  for  images  are  yet  to  be 
discovered,  a few  simple  rules  have  been  found  useful  in  the 
optimization  cf  image  acquisition  and  reproduction  techniques. 

Consider  fcr  example  a black  and  white  reflection  print,  which 
consists  of  a reflective  backing  coated  with  an  eaulsipn  of 
microscopic  grains  of  silver.  The  image  is  formed  by  controlling  the 
amount  of  silver  in  the  emulsion  and  thus  varying  the  relative  light 
absorption  of  the  print,  within  a typical  dynamic  range  of  50  to 
TCO:  1.  Such  a photograph  conveys  its  pictorial  information  to  an 
observer  irrespective  of  illumination  variations  over  perhaps  four  to 
five  orders  of  magnitude.  This  rather  surprising  phenomenon  is  Caused 
by  the  ability  of  the  visual  system  to  "adapt"  to  ambient  levels  of 
lighting  and  thus  to  extract  the  reflection  properties  of  objects 
[1,2].  Studies  cf  the  reproduction  characteristics  of  optimal  images 
[3]  indicate  indeed  that  although  absolute  brightness  influences 
perceived  quality,  the  guality  criterion  within  the  physical 
limitations  of  any  given  reproduction  situation  is  greatly  dependent 
upon  its  ability  to  reproduce  relative  brightness  ratios.  This  fact 
is  intuitively  satisfying  noting  that  pixel  brightness  ratics  are  a 
property  of  the  scene  reflectances  that  is  invariant  to  the  absolute 
intensity  cf  a uniform  illumination. 
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The  implications  of  the  above  visual  phenomenon  ace  that  the 
digital  representation  of  light  intensities  sensed  by  a scanning 
device  should  ideally  be  a measure  of  image  brightness  ratios  rather 
than  arbitrary  absolute  intensity  values.  This  is  easily  implemented 
in  practice  by  recording  the  logarithm  of  the  measured  image 
intensities.  Many  commercially  available  scanners  provide  for  such  an 
option,  usually  called  density  (as  opposed  to  transmittande  or 
reflectance)  scanning.  On  the  reproduction  side,  care  has  then  to  be 
taken  to  preserve  the  recorded  brightness  ratios,  a process  that  is 
facilitated  by  the  inherent  characteristics  cf  the  photographic 
process  to  be  discussed  in  the  next  section. 

The  Photographic  Process:  Exposure  of  a black  and  white  emmlsion 
to  light  and  subsequent  development  produces  a light  absorbing  layer 
characterized  by  its  optical  density  D which  is  defined  as  the 
logarithm  of  the  ratio  of  transmitted  to  incident  light.  With  all 
other  parameters  fixed,  the  optical  density  is  ideally  related  to  the 
intensity  of  the  exposing  light  I by  the  function  [4] 

D = y log  fl  t]  (1) 

where  t is  the  duration  of  the  exposure.  This  function,  well  knpwn  in 
photography,  is  the  Hur ter-Drif f ield  or  D-log  E curve.  actual 
photographic  materials  depart  from  this  idealized  law  at  both  ends  of 
their  useful  dycamic  range.  The  factor  describes  the  "cpntrast"  of 
the  emulsion  and  is  positive  for  an  ordinary  negative  material,  and 
negative  for  a reversal  process.  Because  the  unexposed  emulsipn  and 
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its  substrate  are  not  perfectly  transparent,  an  additional  "fog"  level 
Dq  is  incorporated  into  the  above  equation  yielding 


D = DQ  + Y log  [it]  (2) 

The  light  reflected  from  a print  or  transaitted  through  a slide  is 
related  to  the  incident  light  1 by  [4] 

I = IQ  10-D  (3) 

The  reproduced  light  intensity  I*  is  given  by 

II  = iQ  10"D[lt]  (4) 

Hote  that  if  Y = -1,  the  conditions  for  an  optinun  reproduction  as 
discussed  in  the  previous  section  are  net. 

It  is  not  easy  to  meet  the  relationship  of  eg. (4)  with  actual 
inage  processing  equipment.  Fila  is  typically  exposed  by  a CET,  LED 
or  laser  as  a secies  of  discrete  dots  which  partly  overlap;  the 
exposure  nay  not  be  uniform  over  the  area  of  the  iaage,  etc.  It  is 
possible  though  to  correct  for  such  defects  with  a numerical 
pre-distortion  of  the  digital  iaage  data.  A simple  model,  appropriate 
for  the  correcti.cn  of  a CHI  scanner,  is  discussed  next. 


y 

k 


Calibraticn  of  I/O  Devices:  Actual  image  acquisition  and 
reproduction  devices  have  a number  of  inherent  imperfections  which 
distort  the  final  product.  For  example,  the  measurement  of  pixel 
intensity  in  scanners  is  usually  not  perfectly  logarithmic  .(often 
linear);  the  pixel  intensities  displayed  on  television  monitors  are  a 
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power  function  of  tha  image  signals;  the  light  sensitive  or  light 
emitting  surfaces  of  electron  bean  devices  are  not  perfectly 
honogeneous;  optical  systens  nay  introduce  significant  vignetting, 
etc.  A nunber  cf  procedures  have  been  devised  to  cope  with  such 
inperfections  £5,6].  Por  example,  table  look-up  or  polyoonial 

approximations  nay  be  used  to  correct  for  the  average  deviaticns  of 
the  electro-optical  transfer  function  fron  the  desired  behavionr.  A 
note  refined  (and  expensive)  solution  is  to  vary  the  coefficients  of 
the  correction  as  a function  of  the  gecnetric  inage  coordinates. 

A true  assessaent  of  I/O  device  perfornane  and  the  gathering  of 

| 

physical  data  for  the  design  cf  correction  schenes  is  best  done  by  j 

) 

producing  test  patterns  such  as  step  tablets  and  neasuring  the  optical  j 

density  functions  obtained  on  hardcopy  or  transparency. 

To  illustrate  the  above,  a new  software  correction  technigue  for 
CUT  scanners  is  presented.  It  is  of  nediun  complexity,  but 

computationally  very  fast  and  has  given  excellent  results  with  a CBT  : 

1 

scanner.  The  major  sources  of  distortions  in  this  case  are  ! 

■< 

schematized  in  figure  1.  The  CRT  light  emission  I as  a function  of  i 

the  drive  and  bias  voltages  0 and  UQ  respectively  [7],  as  approximated  j 

■I 

by 

I = fu  + U + U.1YCRT  (5)  j 

U 1 i 

1 

\ 

where  Uj  represents  the  cut-off  voltage  of  the  CRT.  Optical 
vignetting  produces  a darkening  towards  the  image  corners  (figure  2), 
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b)  Numerical  pre- distortion  for  recording  correction 


Figure  4.  6-1.  Distortions  in  CRT  film  recording  and  numerical 
pre- distortion  for  correction. 
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1 

(a)  Constant  brightness  values  photographed  with  a \ 


| Polaroid  camera.  The  darkening  of  the  corners 

I is  evidenced  by  the  small  cut-off  pasted  in  the 

I middle  of  the  photograph. 
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(particularly  annoying  if  one  atteapts  to  produce  a mosaigue,  see 
figure  2b.  Assuming  that  the  vignetting  is  the  only  space-variant 
distortion,  a fast  table  lcok-up  algoritha  has  been  iaplenented,  such 
that  each  source  of  distortion  aentioned  above  is  corrected  for  |n  the 
appropriate  order.  TheYc^T  and  D-log  E correction  of  figure  lb  are 
straight  toward  look-up  tables  based  upon  aeasured  data.  perhaps  the 
aost  interesting  pre-distortion  step  is  the  vignetting  correction. 
Assuaing  circular  syaaetry,  a second  order  polynoaial  of  the  fora 


I'  =1  fA+  B(x2+  y2)] 


P) 


has  been  used  to  boost  the  light  intenities  towards  the  iaage  corners 
where  x and  y are  the  iaage  coordinates  referenced  to  the  screen 
center.  The  values  A/2«-8x  are  stored  in  a one  diaensional  array  C and 
the  correction  is  aade  by  looking  up  this  array  twice  given  the  pixel 
line  and  coluan  indicies  x^  and  y-  . The  results  from  this  fast 
correction  technigue  are  shown  in  figure  5.  The  variations  in  density 
across  a unifcra  surface  are  less  than  0.1  density  units,  whereas  the 
uncorrected  iaage  had  corners  darkened  by  as  much  as  0.35  density 
units. 
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4.7  Spectral  Sensitivity  Estimation  of  a Color  Image  Scanner 
Clanton  E.  Mancill  and  Milliam  K.  Pratt 


The  spectral  sensitivity  of  a color  scanner  must  be  deteriined  in 
order  to  calibrate  its  response.  Direct  spectral  measurements  ever 
the  continuum  cf  the  spectral  band  are  often  difficult  to  obtain. 
However,  responsivity  measurements  can  be  made  through  spectrally 
selective  filters  to  estimate  the  continuous  spectral  sensitivity  of 
the  color  scanner. 


Spectral  Radiance  Estimation: 
multispectral  image  restoration  involve 
radiance  function  c(^)  from  a series  of 


Many  tasks  in  color  and 
the  estimation  of  the  spectral 
observations  of  the  form 


a 
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c(X)  s.(X)d  X + n. 


(D 


where  s.  ( X)  is  the  spectral  sensitivity  of  the  spectral  measurement 
filter  for  i=1,2,...,P  observations.  The  term  represents  additive 
noise  or  uncertainty  in  the  measurement.  Discrete  estimation 
techniques  can  be  applied  to  this  problem  solution  <1>.  The  first 
step  is  to  discretize  the  continuous  integral  to  form  the  vector 
equation 

T 

xi  = ±i  c + n.  (2) 

where  _s.  and  c.  are  Q x 1 vectors  of  quadrature  samples  of  s.  (X)  and 
c(X),  respectively.  Then,  the  set  of  P observations  may  he  arranged 
into  the  P x 1 vector 


x = S£  + n (3j 

T 

where  the  vector  ^ occupies  the  i th  row  of  the  natti*  S.  The 
system  of  equations  represented  by  eg.  (3)  is  normally  highly 
underdeterained  if  safficient  quadrature  mesh  points  are  taken  to 
reduce  the  quadrature  error  to  reasonable  bounds. 

An  estimate  c of  the  true  spectral  energy  distribution  c can  be 
obtained  by  the  generalized  inverse  estimate  <2> 

c_=  S"  x = ST(S  ST)_1  x (4) 

Although  the  generalized  inverse  provides  a minimum  mean  square  error, 
minimum  norm  estimate  of  c,  ill-ccnditioning  of  S coupled  with 


j 
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observational  ducts  can  lead  to  oscillatory  estiaates.  Since  c is 
generally  guite  smooth,  it  is  reasonable  to  impose  some  smopthing 
constraints  on  the  solution,  A common  type  of  smoothing  estimate  is 
given  by  <3> 


c = M'1ST(SM’1ST)'1  x 


(5) 


where  « is  a smoothing  matrix  of  the  typical  form 


M = 


1 -2  1 0 0 0 0 
-2  5-41000 

1-46-4100 
0 1-46-410 

0 01-46-41 


1-46-41  0 

01-46-4  1 

0 0 1 -4  5 -2 
0 0 0 1 -2  1 


(6) 


A third  alternative  is  to  apply  Wiener  estimation  methods  <tt>.  With 
Wiener  estimation,  the  vector  g_  to  be  estimated  is  assumed  to  be  a 
sample  or  a vector  random  process  with  known  mean  and  covariance 
matrix  J?c.  The  Wiener  estimate  is  given  by 


c = K ST(SK  ST+K  ) " 1 x (7) 

— — c— c—  — n — 

where  K is  the  covariance  matrix  of  the  additive  observational  noise 
—n 

assumed  independent  of  c..  As  a convenient  approximation  the 
covariance  matrix  can  be  modelled  as  a first  order  Harkov  process 
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covariance  matrix  of  the  fern 


K = 

— c 


c 

Q 


,Q"1 

0-2 


(9) 


It 

D 0 
Dip 

*Q-1 


vhtre  0 < p < 1 is  the  adjacent  element  correlation  faetpr  and 

represents  the  energy  of  c.  Observation  noise  is  conaonly  modelled  as 

a white  noise  process  with  covariance  equal  to 

2 


n 


K = -jsr—  I 
— n Q — 


<9) 


where  on  is  the  npise  energy  and  I_  is  an  identity  matrix. 


Color  Image  Scanner  Calibration:  A common  problem  in  the 
evaluation  and  calibration  of  color  image  scanners  is  to  determine  the 
total  spectral  response*  of  the  scanner  taking  into  account  the 
spectral  radiance  of  the  illumination  source,  spectral  absorption  and 
scattering  of  the  optics,  and  spectral  sensitivity  of  the 
photodetector.  Direct  measurements  are  often  not  feasible.  Deferring 
to  eq.(1),  let  c (X)  be  redefined  to  represent  the  spectral  sensitivity 
response  of  the  scanner  and  s^  ( be  one  of  P spectral  test  functions. 
The  measurement  procedure  then  proceeds  as  follows.  An  optical  filter 
of  known  spectral  characteristics,  such  as  an  absorption  filter  or 
narrowband  interference  filter  is  introduced  into  the  scanner  and  an 
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filters  whose  peak  transmissivities  span  the  spectral  region  of 
interest.  The  aeasurements  form  the  vector  of  observations,  and  an 
estimation  operation  is  then  invoked  to  obtain  an  estimate  of  the 
scanner  spectral  response. 

In  order  to  evaluate  the  estimation  procedure,  a computer 
simulation  experiment  was  performed  in  which  simulated  measurements 
were  taken  of  a Gaussian  shaped  spectral  function  through  simmlated 
absorption  filters.  Pigure  1 contains  a plot  of  the  spectral  shapes 
of  the  filters.  The  simulated  measurements  were  then  utilized  as 
spectral  observations  for  estimation  of  c{X).  Figure  2 illustrates 
the  performance  of  the  three  estimation  methods  for  simulated 
measurements  through  the  filters.  In  these  experiments  the  mean 
square  fit  between  the  actual  spectral  function  and  its  estimate  was 
least  for  the  simulated  interference  filter  measurements  using  a 
Wiener  estimate  with  P = 0.9  and  a si gnal- to-noise  ratio  of  1000. 

The  spectral  estimaticn  procedures  have  also  been  applied  to  the 
estimation  of  the  spectral  response  of  an  Optronics  Model  S 2000  flat 
bed  scanning  ■ iccodensitometer.  Pigure  3 shows  the  estimate  obtained 
with  absorption  and  interference  filters  for  the  three  estimation 
methods.  No  direct  measurements  are  available  for  the  scanner  so  that 
no  "ground  troth"  can  be  established.  But,  on  the  basis  of  the 
simulation  experiments,  it  is  concluded  that  the  Wiener  estimate 
obtained  with  the  set  of  interference  filters  is  a reasonable  estimate 
of  the  actual  spectral  response. 
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Figure  4. 7-1.  Spectral  shapes  of  absorption  filters 
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4.8  Pedian  Filtering 
W i Ilian  K.  Pratt 

The  aedian  filter  is  a nonlinear  signal  processing  technique 
developed  by  Tukey  <1>  which  is  useful  for  noise  suppression  in 
iaages.  In  one  dimensional  fora,  the  aedian  filter  consists  of  a 
sliding  window  encompassing  an  odd  number  of  pixels.  The  center  pixel 
in  the  window  is  replaced  by  the  aedian  of  the  window  pixels.  The 
aedian  of  a discrete  sequence  a ,a  ,...,a  , for  N odd  is  that  aember 
of  the  sequence  for  which  (N-1)/2  elements  are  smaller  or  equal  in 
value,  and  |N-1)/2  elements  are  larger  or  equal  in  value.  For 
example,  if  the  values  of  the  pixels  within  a window  are 
80,90,200,110,120,  the  center  pixel  would  be  replaced  by  the  valae  110 
which  is  the  aedian  value  cf  the  sorted  sequence  80,90,110,120,200. 


Jn  this  example,  if  the  value  200  was  a noise  spike  in  a moootcnically 
increasing  sequence,  the  median  filter  would  result  in  cpnsiderable 
improvement.  On  the  other  hand,  the  value  200  might  represent  a valid 
signal  pulse  fox  a wide  bandwidth  sensor,  and  the  resultant  image 
would  suffer  seme  less  of  resolution.  Thus,  in  some  cases  the  median 
filter  will  provide  noise  suppression,  and  in  other  cases  it  will 
cause  signal  suppression. 

Figure  1 illustrates  seme  examples  of  the  operation  of  a median 
filter  and  a mean  (smoothing)  filter  for  a discrete  step  function, 
ramp  function,  pulse  function,  and  triangle  function  with  a window  of 
five  pixels.  It  is  seen  from  these  examples  that  the  median  filter 
has  the  usually  desirable  property  of  not  affecting  step  functions  or 
ramp  functions.  Pulse  functions  whose  periods  are  less  than  one-half 
the  window  width  are  suppressed.  Also,  the  peak  of  the  triangle 
function  is  flattened. 

Operation  of  the  median  filtered  can  be  analyzed  to  a limited 
extent.  It  can  be  shown  that  the  median  of  the  product  of  a constant 
K and  a sequence  f(j)  is 

med{Kf(j)]  = K med  (f(j)]  (1) 

Furthermore, 

med  [K  + f(j)  } = K + med  [ f(j) } (2) 

However,  for  two  arbitrary  sequences  f (j)  and  g(j)  it  does  not  follow 
that  the  median  of  the  sum  of  the  sequences  is  equal  to  the  sum  of 
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Figrre  4.8-1.  Examples  of  median  filtering  on  primitive 
signals  - L = S. 
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thtiir  B^dians.  That  is,  it  general 

med  { f(j)  + g(j)l  i med  [f  (j)]  + med{g(j)]  (3) 

The  sequences  80,90,  100, 1 1 C,120  and  80,90,100,90,90  ate  exaaples  for 
which  the  acditive  linearity  property  does  not  hold. 

There  are  various  strategies  for  application  of  the  median  filter 
for  noise  suppression.  Cne  method  would  be  to  try  a median  filter 

with  a window  of  length  3.  If  there  is  no  significant  signal  less, 
the  window  length  could  te  increased  to  five  for  median  filtering  of 
the  original.  The  process  would  be  terminated  when  the  median  filter 
begins  to  do  more  harm  than  good.  It  is  also  possible  to  perform 
cascaded  median  filtering  cn  a signal  using  fixed  or  variable  length 
window. 

The  concept  of  the  median  filter  can  be  easily  extended  to  two 
dimensions  by  utilizing  a two  dimensional  window  cf  some  desired  shape 
such  as  a rectangle  or  a discrete  approximation  to  a circle.  It  is 
obvious  that  a two  dimensional  L x L median  filter  will  provide  a 
greater  degree  of  r.cise  suppression  than  sequential  horizontal  and 
vertical  processing  with  L x 1 median  filters.  But,  two  dimensional 
pcccessing  also  results  in  greater  signal  suppression.  Figure  2 
illustrates  the  effect  cf  two  dimensional  median  filtering  of  a 
spatial  pulse  signal  with  a 3 x 3 square  filter  and  a 5 x 5 plus  sign 
shaped  filter.  In  this  example,  the  square  median  has  deleted  the 
corners,  while  the  plus  median  filter  has  not  affected  the  signal 
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Figure  4.8-2.  Example  of  two-dimensional  median  filtering 
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function. 


Figures  3 and  4 contain  examples  of  the  application  of  median 
filtering  for  iaage  noise  suppression.  In  figure  3 impulse  ncise  uas 
added  to  an  isage.  One  dimensional  aadiaa  filtering  of  length  L=5 
removed  most  of  the  noise  iapulses  mith  only  a snail  loss  in 
resolution.  Alaost  all  errors  vere  removed  for  a median  filter  vith 
L=5,  but  edge  distortion  is  noticeable.  In  figure  4 continuous 
Gaussian  noise  uas  added  to  an  image.  Median  filtering  resulting  in  a 
slight  visual  improvement. 

For  inage  enhancement  applications,  the  median  filter  should 
simply  be  considered  as  an  ad  hoc  tool  for  noise  or  interference 
suppression.  It  should  not  be  used  blindly,  but  rather  its 
performance  should  be  monitored  to  determine  if  its  application  is 
beneficial. 

Reference 

1.  J.  W.  Tukey,  Exploratory  Data  Analysis,  Addiscn-Wesley , 1971. 


(a)  Image  with  impulse  noise 
15  error sAine 


(b)  Median  filtering  of  (a) 
with  L=  3 


(c)  Median  filtering  of  (a) 
with  L = 5 


(d)  Median  filtering  of  (a) 
with  L = 7 


Figure  4.  8-3.  Examples  of  one  dimensional  median  filtering  for 
images  corrupted  by  impulse  noise. 
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*>.  Iaage  Data  Extraction  Projects 
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Iaage  data  extraction  activities  include  the  extraction  and 
measurement  of  image  features,  the  detection  cf  objects  mithin 
pictures,  the  spatial  registration  of  images,  and  the  generation  of 
images  from  one  diaeosional  projections.  Another  facet  of  the  effort 
covers  iaage  pre-processing  operations  which  enable  more  efficient 
machine  data  extraction. 
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5.1  textural  Boundary  Analysis 
William  B.  Thompson 

Previous  reports  have  described  the  development  of  a textural 
distance  function  which  accurately  estimates  the  perceived 
dissimilarity  between  two  textural  regions.  The  textural  distance 
function  model  allous  the  incorporation  of  textural  cues  into  many  of 
the  existing  approaches  to  scene  segmentation.  Texture  may  then  be 
used,  along  with  brightness,  color,  and  any  desired  semantic 
processing  in  determining  object  boundaries.  The  utility  of  textural 
boundary  detectipr  will  be  demonstrated  in  an  edge  criented  system. 

Many  authors  have  developed  edge  finding  systems  which  search  for 
major  discontinuities  in  the  brightness  function  of  the  image  [ 1], 
This  is  normally  dene  by  computing  an  estimite  of  the  derivative  or 
gradient  of  the  iaaga  and  then  finding  the  peaks  in  derivative 
function.  Many  functions  have  been  suggested  for  this  purpose.  A 
ccmmcn  and  often  successful  function  is  called  the  modified  Roberts 
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cross  operator  [2]  and  is  defined  as 


R(i,  j)  = j p(i»  j)  - p(i+l,  j+1)  I + |p(i+l,  j)  - p(i,  j+1)  | (1) 

The  Roberts  "gradient"  is  found  by  sunning  brightness  differenes  in 
two  orthogonal  directions.  Many  more  sophisticated  operators  are 
possible.  In  particnlar,  an  operator  which  returns  edge  orientation 
may  be  quite  useful. 

A procedure  has  been  developed  to  search  for  edges  defined  by 
textural  properties  in  a manner  similar  to  the  Roberts  operator.  At 
specified  intervals  in  the  scene  tc  be  processed,  four  image  regions 
arranged  in  a square  were  considered  (see  figure  1) . The  sum  of  the 
estimated  perceived  textural  differences  between  regions  a and  d and 
between  regions  b and  c was  found.  As  with  conventional  gradient 
operations,  it  was  postulated  that  larger  values  of  this  sum 
correspond  id  tc  textural  edges  cunning  approximately  through  the 
intersection  of  the  four  regions.  In  addition,  an  edge  direction  was 
calculated.  Let  d(i,j)  be  the  computed  dissimilarity  measure  between 
two  regions  i and  j (l(i,j)>0  for  any  two  image  regions).  Theu  a 
textural  boundary  operator  at  the  point  in  the  scene  shown  in  figure  1 
may  be  defined  as 


I 


where  ang  = 0 implies  an  edge  with  negative  slope  at  45  degrees  to  the 
x-axis.  Two  angles  are  possible  since  d |a, d) =d  (b,c)  say  correspond  to 
either  a vertical  or  horizontal  edge.  This  ambiguity  is  straight 
forwardly  resolved  by  considering  d(a,c),  d(b,d),  d{a,b),  and  d(c,d). 

In  the  current  system,  an  edge  map  is  first  produced  by  applying 
the  textural  boundary  operator  at  selected  points  in  an  image.  A 
second  edge  map  is  produced  by  smearing  each  point  in  the  first  map 
along  the  direction  of  edge  orientation.  This  is  doue  to  emphasize 
collinear  edges.  Finally,  actual  edge  points  are  isolated  by  locating 
"ridge  points"  in  the  edge  map.  A ridge  point  is  defined  as  an  image 
point  sufficiently  greater  than  its  neighbors  along  some  direction. 
Much  of  the  code  to  process  the  edge  maps  was  adapted  with  little 
modification  from  a system  originally  designed  to  operate  only  on 
intensity  information  [4]. 

While  most  analysis  systems  designed  to  operate  cn  natural 
imagery  will  use  texture  as  only  one  of  a set  of  multiple  cues  to 
determine  image  organization,  some  way  is  needed  to  evaluate  the 
utility  of  the  textural  boundary  operator  on  its  own.  As  a result, 
this  operator  was  appliad  to  pictures  in  which  the  edges  could  be 
described  as  "purely  textural."  These  test  images  were  created  as 
mosaics  of  textural  patterns  taken  from  pictures  of  natural  scenes. 
Each  component  of  the  mosaic  was  normalized  in  the  same  manner  as  the 

p 

patterns  used  in  the  resolution  experiments.  Thus,  it  was  impossible 
to  distinguish  patterns  based  cn  average  brightness  or  contrast 
criteria. 
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Figure  2 shows  a representative  mosaic  pattern.  Note  that  to  a 
human  observer,  there  are  several  quite  prominent  edges.  Thus,  it  is 
clear  that  human  perception  can  identify  boundaries  on  criteria  other 
than  differences  in  average  brightness.  Figure  2a  is  another  mosaic 
pattern.  Figure  3b  indicates  the  different  textural  regions  present 
in  figure  3a.  In, figure  3a,  a very  prominent  boundary  exists  between 
patterns  a and  b.  The  boundary  between  b and  d is  relatively 
noticeable  while  the  edge  between  a and  d is  hardly  detectable. 
Region  c may  be  viewed  at  cne  level  as  a uniform  teitural  region.  On 
another  level,  however,  the  region  may  be  thought  of  as  being  composed 
of  many  smaller  regions  corresponding  to  the  predominantly  light  and 
predominantly  dark  areas  in  the  pattern. 

The  textural  edge  operator  was  applied  to  these  and  several  other 
mosaic  patterns  asing  several  different  sizes  for  the  basic  blodks  in 
the  operator  (i.e.  the  blocks  in  figure  1).  The  original  mpsaics 
were  256  by  256  picture  elements  in  size.  Pigure  4 is  an  edge  map  for 
figure  3a  using  a basic  blcck  size  of  16  by  16  picture  elements.  No 
post-processing  other  than  the  oriented  smearing  (e.g.  edge  linking, 
noise  cleaning,  etc.)  was  applied.  An  effective  job  has  been  dooe  at 
identifying  the  visually  prominent  boundaries  in  the  mpsaic*  The 
textural  resolution  experiments  would  indicate,  however,  that  it 
should  be  possible  to  achieve  higher  resolution.  Thu3,  it  is  possible 
to  use  olock  sizes  as  small  as  6 or  8 pixels  on  a side.  Figure  2b  is 
an  edge  map  for  figare  2 using  an  8 by  8 basic  block,  size,  all  pf  the 
perceived  boundaries  have  teen  well  located.  Figure  4a  is  an  edge  map 
for  the  mosaic  in  figure  3a  using  the  same  8 by  8 basic  block  size. 


-128- 


(a)  Textural  mosaic  #1 


P 

* 


(b)  Edge  map  for  (a)  using  8x8  regions 


Figure  5.1-2.  Examples  of  textural  mosaics  with 
edge  map. 
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(a)  Edge  map  for  figure  -3a  using  8x8  regions. 


fb)  Edge  map  for  figure  -3a  using  16  x 16  regions. 

Figure  5.1-4.  Edge  map  differentiation  using  8 and  16 
block  regions 
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Again,  the  boundaries  are  well  identified.  The  operator  completely 
degenerates  in  region  c,  however.  A look  at  the  original  picture  will 
show  that  many  of  the  elementary  light  and  dark  areas  are  of 
ccsparable  size  to  the  3 by  8 basic  block.  Thus,  at  this  resolution, 
the  micro- idges  are  a dominant  effect.  This  is  another  example  of  the 
importance  of  realizing  that  perceived  edges  have  a "size"  associated 
with  them  that  is  a function  of  the  size  of  the  objects  being  searched 
for.  Comparable  results  were  obtained  on  the  other  losaic  test 
patterns. 

A difficulty  with  many  of  the  problems  in  automated  image 
description  is  that  it  is  often  almost  impossible  to  quantify  the 
success  of  any  given  approach.  For  example,  the  utility  of  a 
particular  object  isolation  procedure  is  really  only  meaningful  jn  the 
context  of  the  processing  to  follow.  Unfortunately,  the  nature  of  the 
problems  are  so  complex  as  to  make  development  of  completed  systems 
most  difficult.  As  much  of  automated  scene  analysis  involves  the 
simulation  of  perceptual  effects,  the  levelopment  of  lower  level 
operators  described  in  this  report  has  used  human  visual  perception  as 
a performance  gcal. 

The  existence  of  readily  perceived  textural  edges  should  be 
apparent.  In  many  cases,  existing  automated  systems  which  depend  on 
identifying  brightness  discontinuities  will  fail  to  find  these  edges. 
This  report  has  demonstrated  a way  in  which  measures  cf  teytural 
dissimilarity  may  be  incorporated  into  scene  segmentation  systems.  A 
textural  edge  operator  is  developed  which  is  able  to  accurately  locate 
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boundaries  of  a purely  textural  nature. 

The  size  of  the  regiog  over  which  a textural  pattern  is  measured 
has  a significant  effect  on  how  well  that  texture  can  be 
characterized.  Experimental  results  show  that  a dominant  influence  on 
human  textural  resplution  is  the  nature  of  the  patterns  surromnding 
the  region  of  interest.  There  is  a well  defined  trade  off  between 
spatial  resolution  of  a textural  boundary  and  the  ability  to 
distinguish  between  visually  siailar  textures.  The  structural 
interpretation  of  textural  patterns  suggests  several  additional 
methods  for  estimatiag  minimal  resolution  regions.  Unfortunately,  at 
least  one  of  these  measures  (an  auto-correlation  ratio)  is  not 
supported  experimentally.  The  performance  of  the  textural  edge 
operator  for  varying  region  sizes  corresponds  closely  to  the  predicted 
visual  response  frcm  the  rescluticn  experiments. 
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5.2  Image  Segmentation  by  Eomndary  Determination 
Ram  Nevatia 

Finding  boundaries  of  objects  in  an  image  is  a major  concern  of 
scene  analysis.  The  boundaries  ccnsititute  a segmentation  pf  the 
scene.  Conversely,  the  boundaries  may  be  derived  from  a given 
segmentation.  A number  of  segmentation  technigues  have  been  suggested 
in  the  past,  differing  in  their  assumptions  about  the  contents  of  the 
scene  and  in  their  control  structure. 

Using  detailed  specific  knowledge  of  the  objects  likely  fco  be 
present  in  an  image  sixplifies  the  segmentation  process  [1-2],  but 
these  technigues  suffer  from  loss  of  generality.  Another  distinction 
between  various  technigues  is  in  their  control  structures,  such  as 
"tcp-down"  vs.  "bottom-up.”  The  former  treat  an  entire  image  as  one 
object  and  successively  sub-divide  it  into  more  parts  as  needed  [3-4]; 
the  latter  start  frcm  small  atomic  regions  (as  small  as  a single 
pixel)  or  local  edges  and  build  larger  parts  from  them. 

The  bottcm-mp  technigues  are  usually  referred  to  as  being  fedge" 
oriented  or  based  on  "region  growing."  The  edge  based  technigues 
depend  on  detecting  a discontinuity  between  some  properties,  such  as 
brightness  or  color,  of  parts  of  an  image  and  connecting  these 
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discontinuities  to  fora  boundaries.  Region  growing  proceeds  by 
clustering  iiage  points  of  "similar"  properties  in  regions  and  further 
aerging  of  regions  of  siailar  properties  until  a satisfactory 
segmentation  has  been  obtained.  Knowledge  of  isage  properties  has 
been  used  to  guide  the  aerging  of  regions  [5]. 

The  edge  based  approaches  were  initially  used  for  analysis  cf  the 
scenes  of  polyhedral  objects,  the  so-called  "blocks  world."  The 
individual  objects  were  of  unifora,  hcaogeneous  surfaces  and  were  seen 
against  a unifcraly  light  or  dark  background.  Here,  the  edges 
detected  by  a local  edge  operator  usually  correspond  to  the  desired 
object  edges  only.  However,  for  sore  complex  scenes,  the  local 
discontinuities  d.c  not  necessarily  correspond  to  the  object  boundaries 
only;  shadows,  surface  imperfections  and  texture,  and  noise  in  the 
iaaging  devices  being  some  cf  the  causes. 

Consider  the  picture  in  figure  la  showing  a toy  tank  against  a 
background  of  grass.  Note  the  wheels  of  the  tank  are  not  visible  in 
figure  la  because  of  display  linitations.  Figure  1b  shows  the 
intensity  edges  detected  frcn  figure  la,  by  the  application  of  a local 
edge  detector,  known  as  a Hueckel  edge  operator  [6],  at  every  second 
pixel  in  every  ether  rew  of  the  iaage.  This  operator  detects  the 
presence  of  an  edge  in  a circular  neighborhood  and  returns  the 
position  as  well  as  a direction  for  the  edge.  Figure  1b  contains  a 
large  number  of  edges,  most  of  which  dc  not  belong  to  the  desired 
boundary  of  the  tank.  However,  hunans  presented  with  this  edge 
picture  have  no  difficulty  in  perceiving  the  tank.  The  edges  along 
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(b)  Edges  detected  in  (a) 


Figure  5.  2-1.  Edge  detection  for  a picture  of  a toy  tank 


r 


s 
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the  tank  boundaxy  connect  in  a coherent  Nay,  whereas  the  edges  in  the 
grass  region  are  seen  as  being  randoaly  distributed. 

An  algorithm  to  find  groups  of  edges  that  connect  in  an 
approxiaate  straight  line,  to  be  described  later,  is  very  successful 
in  separating  the  tank  boundary  froa  the  background  for  the  above 
exaaple.  This  aethod  of  segaentation  has  the  advantage  cf  being 
general,  as  no  specific  objects  in  the  scene  are  assuaed.  Also,  the 
scheaes  using  texture  properties  defined  over  a regicn  are  sensitive 
to  the  choice  of  the  regicn  size,  and  it  is  difficult  to  locate  the 
boundary  accurately  vithin  a region. 

The  choice  cf  linking  edges  into  straight  lines  was  based  cn  the 
computational  efficiency  of  this  process.  Many  nan-made  and  natural 
objects  have  boundaries  with  elongated  segments.  Further,  any  curve 
can  be  represented  by  piecewise  linear  segments;  the  linking  algorithm 
only  imposes  a constraint  cn  the  aaxiaum  curvature  of  the  segments 
linked. 

Linking  Algccithm:  fluch  work  in  the  past  has  been  concerned  with 
linking  local  edce  elements  into  straight  line  segments.  Two  broad 
classes  of  techniques  are  based  on  the  use  of  the  Hough  transform 
[7-1C],  or  the  use  of  graph  theoretic  methods  [11-12],  However,  these 
techniques  have  been  used  in  situations  where  the  number  of  edge 
elements  is  small  and  most  of  these  elements  belong  to  the  desired 
boundaries.  Their  effectiveness  for  the  problems  considered  here  is 
unclear,  and  in  some  cases  the  computational  costs  are  likely  to  be 
unacceptable  (e.g.  the  algorithms  using  minimal  spanning  trees. 
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require  computational  times  proportional  to  the  cube  of  the  number  of 
edge  elements  to  be  linked).  A detailed  review  may  be  found  in  [13]. 
A description  of  tie  algorithm  developed  follows. 


Pot  this  discussion,  each  edge  element,  ej , is  considered  to  have 
a position  p.  and  an  associated  direction  a.  . Two  oppositely  directed 
edge  elements  are  considered  to  have  different  directions  (differing 
by  180  degrees).  Length  of  an  edge  element,  determined  by  the  size  of 
the  local  edge  operator,  is  unimportant. 


The  entire  360  degrees  range  of  directions  is  divided  in  a number 
of  eguiangular  intervals  (say  12).  Linking  of  edge  elements  along 
directions  in  each  interval  is  examined.  Linking  in  a chosen  interval 
is  constrained  to  edge  elements  having  directions  approximately  within 
this  interval.  The  following  are  the  steps,  in  detail,  for  linking  in 
an  interval  whose  median  angle  is,  say  0.. 


1.  Examine  each  edge  element  and  put  in  a set  E ^ if  the 
lirection  cf  the  edge,  a . is  within  a fixed,  chosen  range,  A0  of 
the  direction  9.  . Note  that  A0  need  not  be  the  same  as  the  width 


of  the  angular  interval.  Figure  2 a shows  the  edge  elements  for 
the  tank  ftjcm  figure  1b,  which  are  within  a 60  degree  range  of 
horizontal  direction  (0^  = 0 degrees). 


2.  Transform  the  co-ordinates  so  that  the  new  x-axis,  lies  along 

0..  Let  (x.  * , y . ' ) be  the  transformed  co-ordinates  gf  the  i-th 
J 1 i 

edge  element  in  set  Ej . 
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3.  Divide  the  image  flane  in  parallel  strips  (buckets)  of  a 
fixed  size  (say  3 pixels  wide),  noraal  to  X*  (figwre  3 shows 
scheaat ical ly  spme  buckets,  with  the  rotated  X-axis  displayed 
horizontally).  Each  edge  element  ei  in  E will  fall  into  one  of 
the  buckets,  determined  by  the  co-ordinate  x^ • . store  the  edge 
elements  in  each  bucket  in  a list  ordered  by  the  value  of  the  y* 
co-ordinate . 


4.  Link  edges  in  each  bucket:  If  two  consecutive  edge  elements 
in  the  edge  list  for  a bucket  differ  in  their  y’  co-ordinates  by 
a distance  ssaller  than  a threshold  TY,  say  2 pixels,  then  the 
two  elesents  belong  to  a common  segment.  e.g.,  tucket  2 in 
figure  3 is  divided  into  segaents  , S^,  S^. 

5.  Link  segments  in  neighboring  tuckets:  If  the  end  points  of 
two  segaents  in  adjacent  buckets  are  within  a distance  cf  T Y in 
their  y*  co-ordnates  and  also  within  a distance  of  TX  in  their  x' 
co-ordinates,  then  the  two  segments  are  merged  into  cne.  Also, 
the  merging  must  not  result  in  a change  of  orientation  of  the 
segment,  e.g.  in  figure  3,  and  3 7,  or  Sg  and  Sg  3re  merged 
but  not  3 ^ and  Sg. 

6.  Retain  only  segments  of  a length  exceeding  a fixed  number 
(say  7) . 


Figure  2b 
elements  of 
lescription  of 


shows  the  linked  segments  resulting  from  the  edge 
figure  2a,  using  the  thresholds  indicated  in  the 
the  algorithm  above.  Figure  2c  shows  linked  segments 
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from  12  intervals  covering  the  entire  363  degrees  range.  Note  that 
segments  from  different  intervals  are  not  linked  though  they  appear  so 
in  the  figure,  and  some  edge  elements  are  connected  to  more  than  one 


f 


segment.  Resolution  of  such  overlaps  and  linking  of  intersecting 
inter-interval  segments  is  straight-forward.  j 

The  above  described  algorithm  uses  many  thresholds  at  various  | 

i 

steps.  However,  the  algorithm  is  relatively  robust  to  these  choices 
and  the  programs  work  well  on  widely  different  scenes  without  changing 
these  thresholds.  The  same  program,  without  change  of  thresholds  has 

i been  tried  on  different  images,  including  the  problem  of  rib  detection 

in  a chest  k-ray,  with  encouraging  results.  The  details  of  the  basis 
of  choice  of  threshclds  are  found  in  [13]. 

Computational  Complexity:  The  various  steps  of  this  algorithm 
reguire  the  processing  of  an  edge  element  either  in  isolaticn  or  in 
comparison  with  its  immediate  neighbors  in  an  ordered  list.  Thus  all 

j 

computing  costs  are  linearly  proportional  to  the  number  of  edges  i 

processed,  except  for  the  possible  costs  of  sorting  the  edge  lists  in 

« 

! step  3 above. 

The  number  of  edges  in  any  single  bucket  is  normally  a small 
proportion  of  the  total  number  of  edges.  Taking  advantage  pf  the 
initial  raster  order  of  the  edges,  the  sorting  time  can  be  limited  to 

< increase  only  linearly  with  the  number  of  pixels  in  the  image.  The 

sorting  details  are  not  discussed  here. 

[J 

[*,  For  the  example  of  the  tank,  the  total  time  to  link  in  12 

\r 

* 

* 

4.- 

V 
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directions  was  20  seconds  on  a PDP-KI10  processor.  The  pcograas  are 
written  in  the  SAIL  language.  The  total  number  of  edges  detected  was 
about  5000.  The  naximum  memory  requirements  were  about  50K,  36  bit 
words. 

The  techniques  described  are  limited  to  discovering  elongated 
segments  of  edge  boundaries.  These  segments  have  to  be  cpnnected  to 
form  complete  object  boundaries.  There  is  sufficient  informatipn  to 
connect  these  segments  as  evidenced  by  our  ability,  as  humans,  to  do 
so  (in  figure  1b  for  example)  without  recourse  to  the  original  grey 
level  picture.  The  segments  cannot  be  simply  connected  to  their 
nearest  neighbors;  some  notion  of  preferred  configurations  is 
required.  Two  lpng  parallel  segments  are  often  boundaries  of  opposite 
sides  of  a part  cf  an  object;  e.g.,  see  the  boundaries  of  the  barrel 
of  the  tank  in  figure  1b.  Information  obtained  by  other  fccms  of 
analysis  of  the  image,  such  as  texture  or  color  analysis,  will  aid  in 
the  connection  of  these  segments.  Alternatively,  these  segments  may 
be  used  to  aid  in  such  analysis. 

References 

1.  R . Bajcsy  and  3.  Tavakoli,  "A  Computer  Recognition  of  Bridges, 
Islands,  Rivers  and  Lakes  from  Satellite  Pictures,"  Proceedings  of  the 
Symposium  on  Hachine  Processing  of  Remotely  Sensed  Data,  Purdue 
University,  Cctoter,  1913. 

2.  J.M.  Tecenbaum,  "On  Locating  Objects  by  their  Distinguishing 
Features  in  Nultisensory  Images,"  Computer  Graphics  and  Image 


Processing,  Vcl.  2,  No.  3/4 , December,  1973,  pp.  308-320. 

3.  R.B.  Ohlander,  "Analys|s  of  Natural  Scenes,"  Ph.D.  Ihesis, 
Computer  Science  Department,  Carnegie  Mellon  University,  Pittsburg, 
Pennsylvania . 


4.  s.  Tsuji  and  F.  Tcaita,  "A  Structural  Analyzer  for  a Class  of 
Textures,"  Computer  Graphics  and  Image  Processing,  Vol.  2,  No.  3/4, 
December,  1933,  pp.  216-231. 

5.  J.A.  Feldman  and  Y.  Yakimovsky,  "Decision  Theory  and  Artificial 

Intelligence:  I.  A Semantics-Based  Region  Analyzer,"  Artificial 

Intelligence,  Vcl.  5,  No.  4,  Winter,  1974,  pp.  349-372. 

6.  M. H.  Hueckel,  "A  Local  Visual  Operator  Which  Recognizes  Edges  and 
Lines,"  Journal  of  the  ACM,  Vol.  20,  No.  4,  October,  1973,  pp. 
(34-647. 

7.  A.K.  Griffith,  "Edge  Detection ^_jji__S4rrpJ.e«  See  ne  s Using  A Priori 
Information,"  IEEE  Transactions  on  Computers,  Vol.  C-22,  No.  4, 
April,  1973,  pp.  271-381. 

8.  R.o.  Duda  and  P.E.  Hart,  Pattern  Recognition  and  Scene  Analysis, 
John  Wiley  and  Sens,  New  York,  1973. 

9.  W.A.  Perkins  and  T.O.  Binford,  "A  Corner  Finder  for  Visual 
Feedback,"  Computer  Graphics  and  Image  Processing,  Vol.  2,  Nos.  3/4, 
December,  1973,  pp.  355-376. 

1C.  S.D.  Shapiro,  "Detection  of  Lines  in  Noisy  Pictmres  Using 


r 


»» 


Second  International  Joint  Conference  on  Pattern 


i 


Clustering, 

Recognition,  August  13-15,  1S74,  Copenhagen,  Danmark,  pp.  317-318.  | 

I 

11.  E.V.  Raaer,  "Computer  Edge  Extraction  froa  Photographs  of  Curved  j 

| 

Cfcjects,"  New  York  University  Technical  Report  CRL-3U,  Deceaber,  1973. 

12.  C.T.  Zahn,  "Graph-lheor etical  Methods  for  Detecting  and 

Describing  Gestalt  Clusters,"  IEEE  Transactions  on  Computers,  Vol. 

C-20,  No.  1,  January,  1 S 7 1 , pp.  68-86. 

13.  R.  Nevatia,  "Object  Boundary  Deterainaticn  iu  a Textured 

Environment, " (to  be  presented)  Annual  ACM  Conference,  October,  1975, 
Minneapolis. 

c.3  Color  Edge  Detection 

Ram  Nevatia  and  Nilliaa  D.  Miller 

A digital  icage  may  be  represented  as  a matrix  of  values  of  a 

function  I(x,y),  defined  at  digitized  points  in  the  image.  For  a 

black  and  white  image,  I is  a scalar  valued  function,  corresponding  to 

the  brightness  of  the  image  at  the  digitized  points.  For  a color 

image,  I is  a vector  valued  function  having  three  components,  say  IR, 

I and  I , the  intensity  values  in  the  red,  green  and  blue  color  bands 
G B 

respectively. 

In  a black  and  white  imaga,  an  edge  is  defined  by  a discontinuity  j 
in  the  scalar  valued  function  I(x,y).  An  eige  in  a color  image  may  be 
defined  in  several  ways.  If  a metric  were  defined  on  the  vector  space  j 


spanned  by  I,  edges  could  be  detected  in  the  new  scalar  space.  Note, 

this  is  similar  to  reducing  the  color  iaage  to  an  equivalent  grey 

level  iaage.  Alternatively,  edges  aay  be  detected  in  the  three 

coaponents  I„,  I _ and  I „ of  I independently  and  a single  edge 
R u D 

determined  from  their  combination.  A scheme  fcr  color  edge  detection 
is  developed  in  the  .following. 

First  consider  the  details  of  edge  detection  in  a single  grey 
level  image.  It  is  useful  to  consider  an  edge  as  having  a position 
and  also  a direction  (a  magnitude  reflecting  the  discontinuity  aay 
also  be  included).  A simple  gradient  operation  followed  by 
thresholding  provides  such  edge  output.  An  edge  is  often  limited  to 
belong  to  certain  classes  of  discontinuities,  e.g.  a step-like  cr  a 
line-like  discontinuity.  Consider  step  edges  only.  Edge  detection 
may  then  be  viewed  as  the  best  fit  of  a neighborhood  of  an  image  by  a 
step  functicn,  and  requires  determination  of  the  position,  orientation 
and  the  magnitude  of  the  step.  Decision  of  the  presence  of  an  edge  is 
based  on  the  size  of  the  step  (and  perhaps  the  quality  of  the  fit) . 

It  was  suggested  by  Binford  £1],  that  a color  edge  be  determined 
by  making  best  fits  to  the  three  functions  XRr  XG  and  Ijj  separately, 
but  constraining  the  orientation  of  the  step  to  be  the  same  for  all 
three  components,  and  the  decision  of  the  presence  of  an  edge  based  on 
the  magnitudes  of  the  three  steps. 

A popular  edge  detector  for  black  and  white  image  has  been 
developed  by  Hueckel  [2].  This  operator  determines  the  presence  of  an 
edge  in  a circular  neighborhood  and  provides  the  position,  prientation 


and  the  magnitude  cf  the  edge.  Briefly,  it  proceeds  by  approxiaat ing 
the  circular  neighborhood  by  expansion  in  a finite  nuaber  of  teras  of 
an  orthogonal  series  of  functions.  Hueckel  claias  the  chosen  series 
to  be  optiaal  under  certain  assuaptions.  L at  ai  be  the  coefficients 
of  the  expansion  for  a given  neighborhood  (i  ranges  froa  0 to  7). 

A best  step  {unction  is  fit  to  the  approxiaated  function  neyt.  A 
step  function,  parameter ized  by  a tuple,  is  expanded  in  the  saae 
series  to  yield  coefficients  (tuple).  The  paraaeters  of  the  step 
are  chosen  to  xiniaize  the  function 


N =y^[a.-  s. (tuple)]1 

1T6  1 1 


An  attractive  part  cf  Haeckel's  approach  is  that  analytic 
solutions  to  this  miniaization  problea  can  be  found,  avoiding 
expensive  searches.  In  particular,  the  orientation  of  the  optiaal 
step  can  be  deterxined  independently  of  other  paraaeters. 

To  extend  this  concept  to  a color  elge,  the  function  to  be 
■inixized  may  be  fcrxulated  as 

n2  = n1+ng+nb  (2) 

The  functions  N„,  N_  and  N are  as  iefined  in  eg.  (1)  for  the  three 
K U 13 

components  of  the  iaage  I,  i.e. 
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(3) 


7 


^[ai  - si (tuple)] 
i=o  lR  1 


2 


where  the  subscript  R refers  to  the  red  component  and  similar 
expressions  exist  for  Nq  and  Ng. 

The  minifixation  process  now  requires  determination  of  three 
tuples  of  parameters,  with  the  constraint  that  all  three  have  the  same 
orientation  parameters.  Again,  it  turns  out  that  the  prientation 
parameter  can  te  determined  independent  of  the  other  parameters. 
Further,  once  the  orientation  has  been  determined,  the  parameters  in 
one  tuple  can  he  determined  independently  of  parameters  in  the  other 
tuples. 

The  algebraic  details  of  the  derivation  are  not  presented  here. 
A black  and  white,  Hueckel  edge  operator  program,  coded  in  assembly 
language,  has  teen  in  use  at  USC  since  last  year.  It  is  possible  to 
use  many  parts  of  this  program,  as  they  are,  in  the  development  of  a 
color  edge  operator.  This  new  program  is  now  being  developed  and 
debugged. 

other  interesting  considerations  for  color  edge  detection  are  in 
the  weightings  of  the  steps  obtained  for  the  three  color  components. 
It  is  expected  that  transformations  of  the  R-G-B  space  to  another 
three  dimensional  space,  which  is  claimed  to  be  Euclidean,  based  on 
models  of  human  perception  developed  at  USC  [3],  should  aid  in  this 
task  . 
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5.4  Image  Boundary  Estimation* 

Nasser  E.  Nahi  and  Kohammad  Jahanshahi 

In  visual  perception,  among  the  most  effective  stimulus 

configurations  are  the  "edges"  outlining  objects  within  an  image,  [ 1 ]. 
This  has  motivated  many  researchers  in  the  area  of  automated  image 
processing,  specifically  scene  analysis,  to  develop  various  techniques 
of  edge  detection  and  boundary  estimation.  An  incentive  for  research 
in  sceue  analysis  is  the  study  of  robotics  £2].  The  available 
information  about  the  shapes  and  sizes  of  physical  objects  constitute 
and  total  visual  intelligence  required  by  a robot.  Such  information 
can  be  provided  through  knowledge  of  object  boundaries. 

The  oldest  method  known  for  boundary  determination  is  that  of 
thresholding  [3].  This  method,  along  with  the  later  procedures  of 
*This  research  was  partially  supported  by  National  Science  Foundation 
ENG  75-03423. 
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are  well  known  ko  be  highly  sensitive 


locating  the  maximum  gradients, 
to  the  sources  of  degradation  phenomenon  [4].  Various  refinements  of 
the  above  methods,  which  tc  some  extent  account  for  the  presence  of 
ncise,  have  been  recently  introduced  [5j. 

In  this  report,  a boundary  estiaator  is  introduced  for  a certain 
class  of  noisy  iaages.  The  iaarges  considered  contain  an  object  of 
interest  within  a background.  Defining  the  set  of  points  which 
separate  the  object  and  the  background  as  "object  boundary,"  a 
recursive  estimator  is  designed  to  yield  an  estiaate  of  the  object 
boundary.  Extensions  cf  the  estiaator  to  multi-object  images  are 
discussed.  The  performance  of  the  estimator  is  illustrated  through 
applications  to  a few  images. 

Problem  Statement:  Consider  the  class  of  images  which  can  be 
partitioned  into  two  regions:  background  and  foregrpund.  The 

foreground  is  assumed  to  form  a "horizontally  convex"  object.  Given  a 
ncisy  version  of  such  an  image,  the  aim  is  to  obtain  an  estimate  of 
the  object  boundary. 

Modeling  cf  Images  by  Replacement  Processes:  An  image  whose  grey 
level  values,  denoted  by  a two-diiensional  function  fc(n,m),  are 


unknown  is 

commonly  modeled  by  the 

given  first 

and 

second 

order 

statistics 

of  b(a,n).  Literature 

in  the  area 

of 

digital 

iaage 

restoration  includes  use  of  this  information,  along  with  a set  of 
observations,  to  derive  a set  of  estimates  (often  a minimum  mean 
square  estiaate)  for  b(m,n)  [6,T],  However,  consistent  in  the  results 
has  been  the  presence  cf  blurry  edges.  Intuitively,  it  may  be 


I 


concluded  that  an  image  model  based  solely  on  the  first  two  moments  of 
b(m,n)  might  be  suitable  cor  reconstruction  of  image  grey  level 
values,  but  it  does  not  carry  sufficient  information  to  adequately 
reconstruct  tbe  object  boundary. 

\ model  for  the  image  signal  b (m,n)  which  explicitly  represents 
the  object  boundary  along  with  tbe  background  and  object  internal 
details  is  giver,  by 

b(m,  n)  = y(m,  n)bQ  (m,  n)  + [ 1-  y (m,  n)]b,(m, n),  O) 

where  bQ  and  represent  the  intensity  values  of  the  object  and  the 
background,  respectively,  and  y carries  the  boundary  information  of 
the  object  withir.  the  image.  The  two-dimensional  functions  bQ(m,n) 
and  bb(m,n)  are  assumed  to  be  sample  functions  of  two  statistically 
independent,  wide  sense  stationary  random  processes  whose  first  two 
moments  are  given.  The  mean  values  of  bD  and  bb  are  indicative  pf  the 
object  and  the  background  brightness  similarities,  whereas,  their 
respective  autocorrelation  functions  are  measures  of  the  object  and 
the  background  textural  information. 

The  binary  valued  function  y(m,n),  another  random  process,  takes 
values  of  1 or  C corresponding  to  the  points  in  the  image  belonging  to 
the  object  or  the  background,  respectively.  In  the  literature,  this 
function  is  usually  known  as  the  image  "characteristic  function1*  £8], 
The  statistical  pcoperties  of  y will  be  described  shortly. 


The  image  acdel  aathe latical ly  represented  by  eg.  (1),  is  tased  on 
a concept  called  a "replacement  process"  where,  by  definition,  a 
segment  of  a function  or  a random  process  is  replaced  by  another 
function  or  random  process  according  to  a certain  rule  [9], 
Considering  that  for  typical  images  the  object  signal,  in  fact, 
"replaces"  a portion  of  the  background  signal,  the  structure  of  this 
model  is  justified.  In  the  model  of  eg.  (1),  replacement  of  the  object 
process  bQ  with  the  background  process  takes  place  according  to  the 
values  of  y . 

For  future  reference,  note  that  the  domains  of  the  sample 

functions  b (m,n)  and  b,  (m,n)  are  defined  to  be  the  entire  image, 
o b 

This  is,  in  fact,  the  main  motivation  behind  introducing  the  concept 
of  replacement  processes  in  the  image  modeling. 

A sequence  of  observations  constructed  as 

y(m,  n)  = b(m,  n)  + v(m,  n)  (2) 

are  assumed  available  for  neasurement,  where  b (a,n)  is  as  defined  by 
eg.  (1)  t and  v(m,n)  denotes  an  uncorrelated  process  representing  the 
observation  noise. 

An  image  scanner  will  now  be  considered  which  transforms  the 
two-dimensicna 1 data  representing  the  noisy  image,  y (m,n) , into 
one-dimensional  data.  The  scanner  output,  in  the  absence  of 
observation  noise,  is  denoted  by  s (k) , where 

s(k)  = Mk)so(k)  + E l-X(k)]  sfa(k)  (3) 
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models  the  image  in  terms  cf  its  grey  level  values  and  object  boundary 
as  viewed  by  the  output  of  a line  by  line  scanner. 

The  structure  of  the  one-dimensicnal  model  of  eq.  (3)  preserves 

the  replacement  processing  concept.  The  functions  sq  (k)  and  (k)  are 

associated  with  b (m,n)  and  b,  (m,n),  respectively.  That  is,  s (k)  and 

o b o 

s^lk)  denote  the  grey  level  values  of  the  scanned  pbjecta  and 

background,  and  are  assumed  to  be  sample  functions  cf  two 

statistically  independent,  cyclo-stationary  randcm  processes  [10], 

whose  first  twc  cements  are  obtainable  directly  in  terms  of  the  first 

and  second-order  statistics  of  b (a,n)  and  b^{m,n)  [6],  As  in  the 

case  of  bo(m,n)  and  b^(m,n),  the  domains  of  the  sample  functions  sQ(k) 

and  s,  (k)  are  the  entire  scanned  image, 
b 

The  binary  valued  process  X (k)  is  the  one  dimensional  counterpart 
of  y(i,n).  Its  statistics  will  be  described  below.  Note  that  the 
statistics  of  X completely  define  those  of  y. 

Let  mj  and  m2  indicate  the  first  and  the  last  lines  of  the  object 
as  viewed  by  the  scanner,  anla^  , ^ represent  the  beginning  and  end 
points  of  the  object  cn  line  l,  respectively.  In  general,  c^ , i*2  , a ^ , 
(3  for  ra  are  random. 

■V  1 “ 

The  function  X(k) , appearing  in  eg. (3) , is  now  defined  in  terms 
of  a and 


•k 

f 
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where  u[  ] is  the  unit  step  function,  J denotes  the  number  of  picture 
elements  in  one  line  of  the  imaje,  and  The  statistics  cf  the 
process  A(k)  can  now  be  given  in  terms  of  the  statistics  of  »^»  »^r 
and  the  sequence 


(5) 


Assume  that  S forms  a first-order  Markov  process.  This 
assumption  is  made  for  the  sake  of  computational  simplicity,  and  it 
emphasizes  the  dependence  of  the  object  boundary  points  on  line  ^ upon 
the  points  located  on  the  previous  line,  -t-1.  It  is  further  assumed 
that  the  require!  density  functions  are  given,  and  that 


P(WJ  W<t-l’nVm2)  =p(Wt- 1’™!* 


(6) 


Notice  also  that 


p(  wt 1 Vi’  mi)  = p(ar  *1 1 V 1»  *1-1’  mi) 


n> 


= p{atlVi'pw'mi)*  p(gJa^V  p-t-rmi 


The  two  diirensipnal  observation  sequence  y(m,n)  in  o<g.  (2)  will 
also  te  replaced  fcy  its  scanned  version  given  by 

y(k)  = s(k)  + v(k)  P) 

where  s (k)  is  as  defined  in  eg.  (3),  and  v (k)  is  a zero  mean  Gaussian 

2 

white-noise  process  with  variance  ct  . 


To  locate  the  object  boundary,  estimates  cf  the  first  and  last 


lines  (aj  and  m2)  and  estimates  of  the  starting  and  ending  points  (a^ 

and  (3  ) of  the  object  are  sought.  The  estimation  procedure  developed 

here,  as  will  be  shown,  requires  the  values  of  s (k)  , and  s (k)  , 

o b 

Kk<N,  where  N is  the  total  number  of  pixels  (picture  elements)  In  the 

image.  Since,  in  general,  these  values  are  not  known  (cases  of  4Ba9es 

with  known  grey  levels  are  exceptional),  the  estimates  of  s (k)  and 

o 

s^(k)  will  be  used  in  their  place.  Such  astimates  can,  for  example, 
be  obtained  by  implementation  of  the  results  in  [11]  where  only 
two-dimensional  statistical  information  on  s(k),  or  y (k)  is  used. 
Notice  that  the  concept  of  replacement  processes  assures  the  existence 
of  the  estimates  of  sQ(k)  and  sb(k)  for  all  Kk<N.  Since  the  ait  of 
this  paper  is  estimation  of  the  object  boundary,  it  will  be  assumed 
that  the  values  sQand  sb  (cr  their  estimates)  are  given. 

The  boundary  estimation  problem,  as  evident  from  e<js.  (3)  and  (8)  , 
is  a nonlinear  estimation  problem.  Furthermore,  due  to  the  type  of 
nonl inearities  involved  (such  as  the  binary  nature  of  Mk)),  the 
available  estimators  based  cn  linearization  concepts  (such  as  extended 
Kalman  filters)  do  not  yield  satisfactory  results. 

In  this  work,  a set  of  maximum  a posteriori  (MAP)  estimates  for 
the  unknowns  B2  * a l ' an^  ^ -t  ace  obtained.  It  is  shown  that  the 

MAE  estimates  will  xinimize  the  following  expression 


Numerical  Derivation  cf  Estimates:  Acjuisition  of  a numerical 
solution  for  the  minimization  process  of  eg.  (9)  is  an  integral  part  of 
this  presentation.  Since  a rigorous  solution  of  eg.  (9) , resulting  in 
a set  of  optimal  estimates  for  ,m2#a^,  and  , is  compmta  tipnally 
unacceptable,  approximate  solutions  are  soujht.  Two  approaches,  shown 
later  to  yield  satisfactory  results,  are  described  in  the  following. 

One  approach  is  to  obtain  the  estimates  of  and  P ^ over  the 
range  m <£<m  , with  the  assumption  that  values  of  m.  and  a_  are  given. 

1 C 1 M 

For  example,  values  cf  and  ra2  nay  be  chosen  as  n^  = 1 and  m2=n, 
implying  that  the  object  boundary  points  lie  on  every  line  of  the 
image.  Then,  if  necessary,  one  may  utilize  additional  structural 
properties  of  the  cbject  to  eliminate  those  boundary  point  estimates 
incompatible  with  the  given  structural  information. 

An  alternative  approach  is  to  consider  the  problem  in  two  steps; 
namely,  solve  fcta^  and0^,  mjCfXm^,  for  a selected  set  of  m^  and  i»2; 
then  solve  for  the  estiaates  of  m^  and  m2  by  replacing  the  estimates 

A A 

a and  ft  for  « and  p . A recursive  procedure  will  result  if  these 
two  steps  are  performed  at  each  scan  line  resulting  in  an  algorithm 
which  yields  a set  of  estimates  for  ai»B2*a{/  and  P^  , concurrently. 

The  farmer  approach  is  computationally  more  attractive.  However, 
it  requires  additional  information,  of  a nonstatistical  geometric 
nature,  on  the  object,  beyond  the  given  statistical  information,  to 
completely  specify  the  object  boundary. 


Computation:  Assume  that  the  first  and  the  last  lines  of  the 


object,  i.e.,  and  n^,  ate  given.  Then,  eg. (9)  can  be  reduced  to 


d 

”nt  [T(w^)-  2a2  In  p(w^  | w^,  ^rr^)]  1 


Now,  from  eg.  (6) 


[H  [T(w^-2a2^np(aJa^i^_i*mi) 


1 = ir^ 
->_2 


Furthermore,  since 


-2a  ^npfpja^a^,  P^njj)]  } 


(4-DJ+p^  (t-DJ+a^-l 

T(w^)  = K(k)-  ^ K(k) 

k=(-L-l)J+l  k=(  -t-l)J+l 


where 


K(k)  = KQ(k)  - lyk) 


then  eg.  (11)  can  te  writter  as 


wn  ^ , 

w [1j  C-2CT  tnp(aJVl’PW’Inl) 

l-  rrij 


-2a2tnp(3tlara^_1,P^_1,m1) 


(l-DJ+P. 


(^-DJ+a^-1 


£ «»  -£ 

k=(-f,-l)J+l  k=(i-l)J+l 


K(k)] } , 
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A recursive,  easily  i iplementable  solution  of  eg.  (14)  is  possible 
if  the  density  functions  of  eg. (13)  are  approximated.  Hence,  the 
Binit izat ion  in  eg. (11)  is  replaced  by 


'{  £ Ch(fi«>  + 8(at,:n 


t = m. 


where 


(-t-l)J+  ot^-1 


g(dt)  "2°  5Z 


k=  (-t-DJ+1 


(*,-i)j+^ 

h(P^)  - -2a2  ^np(^!ara^1J^_1,m1)  + ^ K(k)  (17) 


k=  (<t~l)J+l 


Examples : Several  images  have  been  considered  to  illustrate  the 

results  of  this  section.  Figure  1 depicts  three  such  examples.  All 

the  pictures  have  grid  size  of  256  by  256.  In  each  case  the  msao  and 

variance  of  the  pictures  are  determined,  and  then  a white  Gaussian 

noise  of  specified  variance  is  added  to  each  picture  (figure  2\ 

An  arbitrary  segmentation  procedure  was  performed  to  produde  the 

background,  s,  (k)  , and  foreground,  s (k)  , Kk<256*256,  sample 

b o 

functions  for  each  picture.  The  segsentation  procedure  was  based  on 
replacing  the  object  intensity  values  by  the  maximum  background 


brightness  value  (forming  the  background  sample)  and  the  background 


(a)  Noisy  Square  (J>)  Noisy  Square 

(S/N  = 1.0)  (S/N  = .6) 


(c)  Noisy  Diamond  (d)  Noisy  Diamond 

(S/N  = 1.0)  (S/N  = .6) 
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intensity  values  by  the  minimum  object  brightness  value  (forming  the 
object  sample).  In  general  an  estimator  is  needed  to  perform  the 
segaentation ; however,  since  the  original  images  were  available  here 
(not  usually  the  case),  the  above  technigue  was  a more  convenient 
procedure.  with  values  o£  n ^ and  *2  9i*e*  as  1 and  256,  respectively, 
the  outputs  of  the  boundary  estiaator  are  shown  in  figure  3 
The  signal  to  noise  ratio  (S/N=signal  variance/noise  variance)  of  the 
observed  image  and  the  conjectured  values  of  the  object  maximum  width, 
L,  are  indicated  in  each  figure. 

5.5  Principal  Components  and  Ratioing  for  Multispectral  Image  Analysis 
Guner  s.  Rofcinscn  and  Werner  Frei 

Manual  or  machine  classification  of  multi-spectral  images  is,  in 
general,  made  difficult  by  the  dimensionality  of  the  problem  and  by 
the  fact  that  the  information  of  interest  may  reside  in  subtle 
differences  between  the  spectral  bands.  However,  the  redundancy 
between  multispectral  images  provides  potentiality  for  a reduction  in 
dimensionality  without  an  appreciable  information  loss.  Both  linear 
and  nonlinear  transforms  have  been  studied  to  achieve  such  a reduction 
and  to  enhance  spectral  dissimilarities  for  terrain  classification  of 
the  four  spectral  bands  of  Earth  Resources  Satellite  (ERTS)  imagery. 

The  principal  component  transformation  is  a well-known  linear 
method  by  which  a linearly  independent  (uncorrelated)  set.  of  images  is 
obtained.  The  energy  compaction  property  of  this  transformation  makes 
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L = 100 


(b)  Square  Boundary 
S/N  = 0.  6 


(a)  Square  Boundary 
S/N  = 1.  0 


(c)  Diamond  Boundary 
S/N  = 1.  0 


L = 140 


(d)  Diamond  Boundary 
S/N  = 0.6 


(e)  Girl  Boundary 
S/N  = 10 


(f)  Girl  Boundary 
L =250  S/N  = 0.  9 


Figure  5.4-3.  Boundary  estimates 


it  particularly  attractive  for  the  reduction  of  dimensionality , but 
the  ccmputat ional  load  may  be  considered  excessive  in  some  cases. 

Another  popular  technique  is  to  generate  ratio  images  in  which 
each  pixel  value  is  egual  to  the  rescaled  ratio  of  the  amplitudes  of 
two  spectral  bands.  The  advantage  of  this  ncnlinear  transformation  is 
that  ratios  are  invariant  to  illumination  variations  and 

coxputationally  fast.  The  disadvantage  is  that  there  are  six  possible 
ratio  images  (disregarding  inverses)  with  rather  similar  energy 
contents. 

Principal  Component  Analysis  of  Hultispectral  Images:  Principal 

components  analysis  of  ERTS  bands  is  motivated  by  the  desire  to 
extract  the  most  significant  spectral  components  from  the  available 
four.  This  dimensionality  reduction  also  results  in  preserving  most 
of  the  ERTS  information  in  a smaller  number  of  components. 

The  principal  component  analysis  of  ERTS  data  involves  finding  a 
unitary  transformation  matrix  which,  when  applied  to  the  four  bands, 
results  in  a new  set  of  bands  (principal  components)  having  several 
desirable  characteristics:  the  principal  components  are  uncorrelated 

and  each  component  has  a variance  less  than  the  previous  component. 

The  principal  components  are  obtained  from  the  original  four 

matrix  multiplication 

y = A x (1) 

of  spectral  intensities  on  four  ERTS  bands*  y is 


spectral  bands  by  the 


where  x is  the  vector 


the  vector  of  principal  components  and  A is  the  4x4  Karhunen-*Loeve 

transformation  matrix.  This  matrix  is  derived  by  diagonalizing  the 

spectral  covariance  matrix  c of  the  spectral  bands.  The  rows  of  A 

— x 

are  the  normalized  eigenvectors  of  C . The  covariance  matrix  of  the 

— x 

principal  components  is  then 


T 

C = A C A 

-y  — x - 


Ml)  0 0 0 

0 X(2)  0 0 

0 0 X(3)  0 

000  X(3) 


(2) 


where  x^  , X2  , X^  , and  (the  variances  of  the  principal  components) 
are  the  eigenvalues  of  ordered  such  that  X^  >X^>X3>X^. 


It  should  be  notad  that,  since  i is  a unitary  transformaticn,  the 
total  data  energy  is  invariant.  That  is 


i = l i = 1 


(3) 


where  the  cr.  , are  the  variances  of  the  original  ERTS  bands.  As  an 
example,  figure  1 shows  four  ER1S  images,  and  figure  2 presents  the 
principal  ccmpcnents  planes.  All  images  have  been  enhanced  by 
histogram  manipulation  before  display.  The  spectral  covariance  matrix 
C of  the  four  ERTS  bands  is  obtained  by  computing  the  spectral 
covariance  matrix  on  64  x 64  blocks  of  EHTS  pictures,  (each  512  x 512 
pixels)  and  then  averaging  over  all  the  blocks.  Exhibit  1 contains 
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Exhibit  5.  5-1 


< 


Statistical  Data  on  Principal  Components  of  ERTS  Planes 


spectral  covariance  matrix 


57.  16  75.80  39.23 
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the  measured  ERTS  covariance  matrix,  the  computed  covariance  matrix  of 


the  principal  components  planes,  and  the  corresponding  eigenvalues. 
It  should  be  noted  that  the  first  two  principal  components  represent 
971  of  the  total  energy. 

Band  Ratios:  Ratioing  of  ERTS  pictures  is  a useful  pre-processing 
technique  for  multispectral  recognition  and  classification. 

Signatures  obtained  from  a training  sample  under  one  set  of  conditions 
may  not  have  a good  discrin inaticn  capability  for  a given 

classification  scheme  if  the  same  area  is  observed  un  ler  a different 
set  of  conditions.  If  the  changes  result  from  simple  multiplicative 
factors  such  as  the  brightress  level,  then  the  ratio  of  the  bands  will 
be  invariant. 

Tailing  various  ratios  of  the  green,  red  and  the  two  infrared 
hands  (bands  4,  5,  6,  and  7,  respectively)  of  the  ERTS  data  results  in 
elimination  of  brightness  variations  due  to  topographic  relief.  Such 
ratio  images  have  been  shown  to  be  more  useful  for  determining 
boundaries  between  lithologic  units  and  vegetation  groups  £1],  Ratios 
may  be  taken  to  emphasize  variations  due  to  color  also.  Such  raticing 
processes  produce  a color  display  whose  color  variations  are  more 
indicative  of  material  variations  than  the  simple  pseudocolor 
displays. 

Ordinarily,  ratio  images  are  obtained  by  formirg  a scaled  ratio 
two  hands,  (direct  ratio).  Logarithmic  ratio  images  are  produced 
tpplyinj  a logarithmic  stretch  to  a ratio  image.  The  advantage  of 
1 jjirithaic  ratio  is  a greater  tolecanca  to  quantization  error. 
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In  the  cases  studied,  it  has  appeared  that  logarithmic  ratio 
images  contain  more  visual  information  than  direct  ratio  images.  It 
is  felt  that  experiments  with  more  images  are  necessary  to  confirm  the 
above  conclusicn. 

As  an  example,  figure  3 shows  the  logarithmic  ratios  of  the  E RTS 
pictures  shown  in  figure  1.  These  ratio  images  have  been  enhanced 
using  the  same  histogram  manipulation  algorithm  as  the  original 
images.  The  choice  of  ratio  images  for  a certain  classification 
scheme  depends  on  the  data  and  the  application. 

The  covariance  matrix  of  various  ratios  could  give  some  insight 
in  choosing  a set  cf  ratios  for  a classification  scheme:  ratios  that 
are  uncorrelated  are  likely  to  produce  better  results  than  those  that 
are  highly  correlated.  This  idea  suggests  the  use  of  the  principal 
components  of  ratios  instead  of  ratios  themselves.  Exhibit  2 contains 
the  normalized  covariance  matrix  and  eigenvalues  of  the  logarithmic 
ratios.  It  is  observed  that  the  first  two  or  three  principal 

components  contain  most  of  the  relevant  information  in  ratio  images. 

i 

. This  can  also  be  verified  by  studying  the  principal  components  shown 

in  figure  4. 

Reference 
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Processing  to  Regional  Geologic  Problems  and  Geologic  Happing  in 
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Exhibit  5.  5- 2 

Statistical  Data  on  Principal  Components  of  ERTS 
Logarithmic  Ratio  Planes 
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Image  Processing  Systems  Projects 


The  following  describes  the  image  processing  systems  projects 
which  are  concerned  with  the  development  of  image  processing  hardware 
and  software  systems. 

6.1  Hardware  Projects 

Toyone  Mayeda 

A real  time  color  imace  magnetic  tape  recorder/playback  system  is 


under  development.  The  recorder  is  to  be  used  to  record  real  time 
digitized  television  signals  at  a 600  ips  rate  and  played  back  at  a 
1-7/8  ips  rate  to  transfer  the  data  to  the  PDF-10  computer.  The 
inverse  process  is  performed  to  produce  real  time  television  signals 
from  coded  computer  records. 

Delivery  of  the  Emerson  (Orion)  digital  magnetic  tape 
recorder/play tack  unit  has  been  delayed  due  to  difficulty  in  meeting 
the  bit  error  rate  and  maximum  skew  specifications.  Emerson  is 
presently  redesigning  the  tape  transport  mechanism  to  reduce  the 
problem.  It  is  also  planned  to  increase  the  deskew  buffer  capacity  in 
the  interface  hardware  which  was  developed  at  USC.  Delivery  is  now 
planned  for  1 January  1976. 

A second  digital  image  television  display  system,  which  is  being 
developed,  is  presently  in  the  check  out  and  testing  phase.  This  unit 
receives  digital  picture  data  from  the  ARPANET,  acting  as  a virtual 
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TIP  terminal,  and  produces  a modulated  television  signal  Cor 
connection  to  the  antenna  terminals  o£  an/  commercial  television 
receiver. 

6.2  Software  Projects 
Dennis  Smith 

The  software  effort  cf  the  Image  Processing  Institute  ( X PI) 
programming  group  has  been  centered  on  tvc  projects.  The  first  has 
been  the  implementation  of  a network  of  mini-computers,  and  the  second 
the  augmentation  of  the  library  of  image  processing  user  programs. 

The  purposes  of  the  network  of  mini-computers  are  to  handle 
communication  among  the  larger  computers  of  the  Engineering  Computer 
Laboratory  and  the  Image  Processing  Laboratory,  and  between  these 
computers  and  machines  at  other  sites,  and  to  handle  lowest  level 
protocols  with  image  processing  devices. 

The  primacy  advantage  of  this  network  is  the  freeing  of  the 
larger  computers  from  the  task  of  minutely  supervising  complex 
devices,  many  of  which  cause  fcegueut  interrupts  that  are  demanding 
upon  a processor's  time.  All  communications  among  the  larger 
computers,  and  hetween  them  and  the  specialized  devices  are  carried  on 
in  message  packets  which  are  blocks  of  data  that  can  be  passed  about 
with  a minimum  c£  interrupts. 

A second  advantage  of  the  network  is  one  of  reliability.  Should 
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the  PDP-11  mini  which  is  controlling  a key  device  become 
ncn-opard t icna 1 , the  software  for  that  device  can  be  easily  moved  to 
another  mini,  the  device  plugged  into  that  machine,  and  service 
restored.  Should  the  PDP-10,  the  principal  computer  for  user 
software,  be  unavailable,  the  LIP-2100  cr  the  IBM  36C/44  can  be  used  in 
this  capacity,  as  the  user  software  is  written  in  portable  FCRT  RAN. 

To  date,  the  two  programs  which  will  run  on  all  the  11's,  the 
supervisor  program,  and  the  network  control  program  (NCP) , which 
manages  the  routing  of  message  packets  from  the  source  to  the 
destination  computer,  are  both  completed.  Remaining  to  be  finished 
are  the  service  programs  to  handle  each  of  the  image  processing 
devices  on  the  11's,  and  the  NCPs  for  each  of  the  larger  computers. 

The  second  area  of  concentration  is  user  software.  Several 
personal  programs  of  the  IPI  faculty  and  staff  were  obtained  frpm  the 
individuals  who  wrote  them  and  were  added  to  the  IPI  library  after 
modification  to  make  them  more  useful  to  the  general  community.  All 
of  the  following  were  standardized  to  conform  to  parameter  input 
conventions  of  the  other  library  programs,  and  generalized  to  process 
images  which  are  any  power-of-two  size  smaller  than  or  equal  to  1024. 
All  programs  run  in  an  interactive  mode,  asking  the  users  questions  as 
to  what  he  wants  done.  These  programs  are  described  below. 

CONVOL  - a program  fcr  performing  two  dimensional  convolution  was 
generalized  to  provide  a choice  of  impulse  response  arrays  (or  allow 
the  user  to  enter  his  own)  in  sizes  3x3,  5x5,  or  7x7. 
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HSTMOD  - a program  to  perform  notification  of  histograms,  which  j 

J 

does  egua  1 izat icn,  exponentiation,  cr  a "gamma"  function  upon  the 
histogram  of  a picture. 

PICOUT  - a program  for  contrast  manipulation  was  expanded  to 
perform  the  following:  clipping,  labeling,  reformatting  (paaking, 

unpacking,  integer-real,  rea 1- integer) , and  application  of  one  of  a 
variety  of  transfer  functions:  positive  linear,  negative  linear, 

sawtooth,  slicer,  eye,  half  power,  third  power,  log,  or  a user-defined 
step  function  (256  steps)  with  automatic  scaling, 

WEDIAN  - a median  filtering  program  which  offers  three  choices  of 
filtering:  N1S,  which  computes  the  median  for  each  position  of  a 

rectangular  window  as  it  scans  the  picture  file;  NEDX , which  computes 
the  median  for  each  position  of  a cross  window;  and  MOVAVG,  which 
computes  the  mean  fox  each  positon  of  a rectangular  window.  All  of 

the  above  may  be  used  with  any  window  size  1 x 1 to  1 1 x 1 1 . 

CFIL  - a program  to  dc  image  restoration  and  wiener  filtering. 

It  allows  specification  of  a blur,  correlation  coefficient,  and 
signal-to-noise  ratio,  and  an  inplulse  response  matrix  up  to  31  x 31. 

6.3  A Synthesis  Procedure  for  Optical  Honlinearities 
Stephen  R.  Dashiell  and  Alexander  A.  Sawchuk 

A general  technigue  for  implementing  nonlinear  nonmonotonic 
function  incoherent  optical  parallel  signal  processing  systems  has 
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been  .inscribed  in  recent  publications  [1-4J.  The  technique  operates 
by  using  special  halftone  screens  and  high  contrast  (binary)  optical 
input  devices  to  effectively  pulse-width  modulate  the  input.  The 
selection  cf  diffraction  orders  in  a Fourier  transform  produces  a 
desampled  output  which  is  a point  nonlinear  function  of  the  input. 

A very  complete  analysis  of  the  entire  process  has  been  performed 
[4].  One  generalization  that  has  been  found  is  that  the  halftone 
prefiles  (cells)  theaselves  which  determine!  the  dot  size  b need  not 
be  monotonic.  Thus,  the  effective  periodicity  of  the  preprocessing 
can  be  changed.  The  effect  is  to  reduce  the  diffraction  order 
necessary  to  achieve  nonmonotonic  operation.  So  many  design  variables 
are  now  available  that  the  class  of  mathematical  operations  possible 
and  ease  of  icplementaticn  has  been  greatly  extended. 

An  exact  synthesis  procedure  for  nonlinearities  using  ordinary 
mcnotonic  cells  has  been  made  and  is  summarized  here  for  the  case  of 
linear  on  -dimensional  scenes.  emitting  wavelength  and  geoieftrical 
factors  for  clarity,  the  general  expression  for  the  amplitude  in  the 
transform  plane  resulting  from  an  infinite  grating  of  opague  bars  of 
width,  b,  and  period,  a,  with  unit  amplitude  illumination  is 


J{y*>l=6(fx>-X)  si"'  'T1 


where  the  y dimension  is  suppressed.  Ey  selecting  these  diffraction 
orders  with  simple  spatial  filters,  the  sine  terms  in  eg.(1)  indicate 
that  nonmonotonic  behavior  can  be  expected.  In  the  special  case  of  a 


zero  order  <n=0)  selection,  an  intensity  output 


*out(0)  = (1  - bW 


is  expected  free  eq.(1)  after  inverse  Fourier  transfprming  and 
squaring.  Por  a first  order  (n=1)  selection,  the  output  intensity 


r 1.2  nb  . 

1 4.n\  = sin  ) 

out(l)  2 a 

TT 


a function  which  is  ncnacnctonic  in  b. 


Because  of  the  halftone  process,  the  value  of  b in  these 
expressions  is  a function  of  the  continuous  input  intensity  I . A 
one-dimensional  halftone  screen  can  be  described  as  periodic 
symmetrical  cells  centered  on  x-0  and  extending  from  -a/2  to  a/2,  each 
with  a density  function  f (x) . The  intensity  I transmitted  by  the 
input-screen  combination  in  each  cell  is 


I = I.  10 
t in 


and  this  function  is  imaged  or  contact  printed  onto  a binary  clipping 
medium  with  effective  cutoff  I*.  Since  there  is  no  exposure  if  I s:  I • 
and  full  exposure  if  I >!',  opaque  bars  result  where  x is  ssch  that 


I.  a I'  iof(x) 


Taking  logarithm  yields 


1% 

If 


* 


: ** 

\* 


bl 


* 5 f"’  [‘°810  i~v) ] <6’ 

where  f * is  the  inverse  of  f.  The  cell  siza  b is  simply  twice  x for 
halftone  cells  symmetrical  about  tfce  spacing  a. 

Combining  eg.  (6)  with  egs.  (2)  or  (3)  for  the  appropriate  order 
gives  the  overall  mapping 

= «o'V=  11  - 2f'1[logio<^1>]/a>2  171 

for  the  zero  (n=0)  order,  and 

'out  = 8lV  = (-“‘"(f1  ' f_I  [1Og10  7E])/i’) 

for  the  first  order  (n=1).  Similar  expressions  can  be  obtained  for 
two-dimensional  cells  and  various  selections  of  diffraction  orders. 

These  expressions  for  transforms  and  dot  sizes  ate  valid  only  in 
local  regions  of  constant  input  values.  Input  inf ormaticn  produces 
low  spatial  frequency  modulation,  and  the  complete  expression  for  the 
transform  is  much  more  complicated.  The  halftone  process  assumes 
input  samplinj  at  a rate  sufficient  to  avoid  aliasing,  and  these 
results  describe  the  local  nonlinear  effects  if  desampling  filters 
choose  the  lew  frequency  input  information. 
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The  following  procedure  can  be  used  to  obtain  the  cell  profile 
f (x)  and  diffraction  crder  for  one-dimensional  screens  given  a desired 


f 

* 

h 

ir, 

(' 
l 1 

i 


I 


= h (I . 
out  1 in 


) : 


I.  Determine  the  minimum  diffraction  order  n to  be  used  by 
counting  the  number  of  sign  changes,  q,  in  the  slope  of  the  transfer 
function.  If  q is  zero  and  the  initial  slope  of  h is  negative,  the 
n=C  order  can  be  used  directly.  If  q is  zero  with  positive  initial 
slope,  the  n=1  crder  must  be  used.  For  q greater  than  zero,  add  one 
to  g if  the  initial  slope  is  negative  to  obtain  g'.  If  q is  greater 
than  zero  and  the  initial  slope  is  positive,  then  q'=q.  The  number  of 
slope  changes  in  the  general  IQut^nj  versus  b curve  is  given  by  2n-1 
for  h>0,  thus  n is  selected  so  that  q*  is  less  than  or  equal  to  2n-1. 
This  procedure  determines  the  minimum  n,  so  that  a larger  order  can  be 
used  if  desired. 


II.  Normalize  the  desired  function  by  scaling  so  the  largest 

Iout  equals  the  maximum  for  the  particular  crder  used.  For  n = 0, 

I <1;  for  n>C,  I <1/n 
out  out 

III.  Equate  h (Iin  ) with  the  appropriate  general  expression 

g (I.  ) of  the  form  eg (10)  or  eg.  (11)  for  the  particular  order  n 
n in 

used.  Solve  this  equation  for  f llcg^g  (I.  /I')). 

IV.  Solve  for  f(x)  by  selecting  a solution  such  that  f (*)  is 

mcnctonic  and  the  initial  slopes  of  h(I.  ) and  g (I  ) have  the  same 

in  n m 

sign.  Whenever  the  slope  of  h(I.n)  changes  sign,  the  halftone  cell 
size  must  abruptly  increase  so  that  tfce  diffraction  output  remains  the 
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f 


same  while  jumping  tc  a region  of  g (I.  ) c£  opposite  slope. 

n m 

An  example  of  this  procedure  is  the  synthesis  of  an  optical  level 
slices,  or  intensity  bandpass,  with  the  characteristics  shown  in 
figure  1.  This  function  is 


I = h(I.  ) = K,  1*1.  M 
out  in  cl  m cl 


0, 


otherwise 


(9) 


and  it  has  one  sign  change  in  slope,  so  g equals  one.  The  initial 

slope  is  positive,  so  q*  is  one  and  the  first  (n=l)  diffraction  order 

2 

can  be  used.  Normalizing  the  function  h(I.  ),  gives  K equal  to  1/rr  , 

m 

and  equating  with  g ^ (I  in)  gives 


f 1(log10Iin)  = (a/2rr)sin  1 


(n[  h(I.  )] 
m 


(10) 


where  the  clip  level  I''  is  assumed  unity  for  simplicity.  For  I <1 


f_1(log  I.  )=  ( --)  sin_1(n[0]2  ) C11) 

10  m ztt 

= 0 or  a/2 

Selecting  the  zero  solution  to  satisfy  the  aonotonic  cell  condition, 
results  in 

f(0)  = logjo^i  <12> 
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(13) 


For  Icl<Iin<Ic2 


f'1(,ogioW  = 

n 


and  solving  gives 


« t > = '°holc2  ,14» 

as  a point  of  discontinuity  of  f.  For  Ic2<Ijn 

f_1(log  -I.  ) = (~  ) ain”^(rrC  0 ]^)  = 0 or  a/2  (15) 
lu  in  c,  tt 

Here  the  a/2  solution  is  selected  to  satisfy  the  aonotcnic  cell 
condition.  Ihis  is  the  end  point  of  the  profile  having  period  a. 
This  function  f (jc)  , 0<x<a/2  shown  in  figure  2,  has  been  experiaentally 
demonstrated  [2-3].  The  width  of  the  level  is  controlled  by  the  step 
size  in  f(x),  and  the  level  location  is  controlled  by  the 
preprocessing  step.  In  general,  the  halftone  cell  profile  aay  cpabine 
sacoth  and  discontinuous  functions,  leading  to  transfer  functions 
h (Iin  ) with  both  s aooth  and  Uniting  nonlinear  characteristics. 

The  analysis  of  systea  effects  due  to  low  contrast  (finite  gaaaa) 
input  aedia  is  well  underway.  In  the  zero  order  the  aajor  effect  is  a 
change  in  the  transfer  function;  in  the  first  order,  this  effect  is 
coabined  with  a less  of  diffraction  efficiency.  These  effects  are  not 
serious  in  practice,  and  scue  techniques  of  pre-coapensating  halftone 
cells  to  correct  for  low  gaaaa  have  been  developed.  These  appear  very 
proaising  for  practical  iupleaenation,  particularly  with  real-tiae 
input  devices.  A series  of  coaputer  routines  have  been  written  to 


-183- 


iteratively  synthesize  cell  profiles,  produce  input/output  curves  and 
study  effects  of  paraaater  variation  on  the  results. 


A nunber  of  experimental  halftone  screens  have  been  sade  and 
tested.  A computer-controlled  optronics  flatbed  microdensitometer  has 
been  used,  and  direct  plots  on  highly  resolution  film  have  been 
adequate  to  make  gopd  quality  screens.  Host  of  the  screens  have  been 
one-dimensional  line  gratings,  and  plotting  aperture  sizes  down  to 
10  m.  have  been  used,  Kodak  50-427  sheet  film  is  used  for  the  screens 
because  of  its  high  resolution  (>250  lines/mm.)  and  good  line  bolding 
ability.  Seme  of  the  functions  which  have  plotted  and  tested  Mith 
good  results  so  far  include:  intensity  level  slicers,  intensity  notch 
filters,  logarithms,  and  exponentials.  Experimental  verification  of 
other  functions  is  underway. 
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